7Консервативные системы с одной степенью свободы

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

H (x, y)

— какая-то гладкая функция,

H (x, y) : R^{2} \to R

. Система

\begin{matrix} ⎧ ⎨ ⎩ \begin{matrix} ˙ x = \frac{\partial H}{\partial y} ˙ y = - \frac{\partial H}{\partial x} \end{matrix} \\ (7.1) \end{matrix}

называется системой в гамильтоновой форме. Функция

H

называется функцией Гамильтона

Утверждение 1. Функция

H

является первым интегралом системы (7.1)

Доказательство. Посчитаем производную функции

H

вдоль векторного поля

v = (\frac{\partial H}{\partial y}, - \frac{\partial H}{\partial x}) .

Имеем:

L_{v} H = \frac{\partial H}{\partial x} ˙ x + \frac{\partial H}{\partial y} ˙ y = \frac{\partial H}{\partial x} \frac{\partial H}{\partial y} - \frac{\partial H}{\partial y} \frac{\partial H}{\partial x} = 0

Иными словами, скорость изменения функции

H

при движении вдоль решений дифференциального уравнения (7.1) равна нулю, то есть

H

— первый интеграл. (См. утверждение 2 из главы 6).∎

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

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

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

\begin{matrix} ¨ x = F (x) . \\ (7.2) \end{matrix}

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

F

, которая зависит только от положения частицы.

import matplotlib.pyplot as plt
import numpy as np
import qqmbr.odebook as ob
# see https://github.com/ischurov/qqmbr/blob/master/qqmbr/odebook.py

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

Рис. 7.1: На точечное тело с координатой

x

действует сила

F (x)

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

\begin{matrix} {\begin{matrix} ˙ x = y ˙ y = F (x) \end{matrix} \\ (7.3) \end{matrix}

Можно найти фазовые кривые этого уравнения

\begin{matrix} \frac{d y}{d x} & = \frac{F (x)}{y} y d y & = F (x) d x \int y d y & = \int F (x) d x \frac{y^{2}}{2} & - \int F (x) d x = C \\ (7.4) (7.5) (7.6) (7.7) \end{matrix}

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

$U (x) = - \int F (x) d x$ — потенциальная энергия.
$V (y) = \frac{y^{2}}{2}$ — кинетическая энергия.
$H (x, y) = U (x) + V (y)$ — полная энергия.

Утверждение 2. Функция

H

является первым интегралом системы (7.3).

Доказательство. Во-первых, это следует из (7.7). Во-вторых, можно заметить, что система (7.3) — гамильтонова. Наконец, можно явно посчитать производную вдоль векторного поля:

L_{(y, F (x))} H = \frac{\partial H}{\partial x} y + \frac{\partial H}{\partial y} F (x) = - F (x) y + y F (x) = 0

Утверждение доказано.∎

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

Пример 1. Уже знакомое нам уравнение осциллятора.

F (x) = - x .

Его потенциальная энергия:

U (x) = \frac{x^{2}}{2}

, полная энергия

H (x, y) = \frac{x^{2}}{2} + \frac{y^{2}}{2}

. Построим график полной энергии:

import plotly
import plotly.graph_objs as go
from plotly.offline import iplot as plot
from plotly.offline import init_notebook_mode

init_notebook_mode()

import numpy as np

X, Y = np.mgrid[-2:2:100j, -2:2:100j]
Z = X**2/2 + Y**2/2

data = go.Surface(x = X, y = Y, z = Z)
plot([data])

Рис. 7.2: График полной энергии гармонического осциллятора

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

import matplotlib.pyplot as plt
import numpy as np
import qqmbr.odebook as ob
# see https://github.com/ischurov/qqmbr/blob/master/qqmbr/odebook.py

plt.figure(figsize=(6,6))
ob.axes4x4(labels=('x','y'))
F = lambda x, y: y**2 / 2 + x**2 / 2
ob.mcontour(np.linspace(-5,5,300),np.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=plt.cm.coolwarm)
ob.vectorfield(np.arange(-4,4,0.5),np.arange(-4,4,0.5), lambda x, y: (y, -x))

Рис. 7.3: Линии уровня функции полной энергии гармонического осциллятора

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

(0, 0)

Пример 2. «Перевернутый маятник» вблизи положения равновесия.

F (x) = x .

Здесь сила пропорциональна отклонению и действует в сторону увеличения отклонения.

Потенциальная энергия $U (x) = - \frac{x^{2}}{2}$ , полная энергия $H (x, y) = \frac{y^{2}}{2} - \frac{x^{2}}{2}$ . Как и в предыдущем примере, построим график и линии уровня полной энергии.

import plotly
import plotly.graph_objs as go
from plotly.offline import iplot as plot
from plotly.offline import init_notebook_mode

init_notebook_mode()

import numpy as np

X, Y = np.mgrid[-2:2:100j, -2:2:100j]
Z = -X**2/2 + Y**2/2

data = go.Surface(x = X, y = Y, z = Z)
plot([data])

import matplotlib.pyplot as plt
import numpy as np
import qqmbr.odebook as ob
# see https://github.com/ischurov/qqmbr/blob/master/qqmbr/odebook.py

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

Рис. 7.4: График и линии уровня полной энергии перевернутого маятника вблизи положения равновесия

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

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

Определение 2. Математическим маятником называется механическая система, состоящая из стержня, один конец которого закреплён, а к другому прикреплен маленький грузик (материальная точка). Стержень может свободно вращаться вокруг закрепленного конца. Трение отсутствует. Масса стержня считается нулевой.

Рис. 7.5: Силы, действующие на маятник. Автор изображения: Krishnavedala / Wikimedia Commons, лицензия CC BY-SA. Оригинальное изображение.

Пусть маятник отклонён от положения равновесия на угол

θ

. Рассмотрим силы, действующие на грузик (см. рис. 7.5). На него действует сила тяжести

m g

(здесь

g

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

m g cos θ

) и перпендикулярно ему (

m g sin θ

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

m g cos θ

. Интересующая нас проекция силы тяжести направлена перпендикулярно к направлению стержня и равна

m g sin θ

: её никто не компенсирует и именно она приводит маятник в движение.

Теперь нам нужно записать соответствующее дифференциальное уравнение с помощью второго закона Ньютона. В качестве фазовой координаты мы рассматриваем угол $θ$ , но положение (а значит и скорость, и ускорение) грузика зависит не только от $θ$ , но и от длины маятника $l$ . Нетрудно видеть, что скорость и ускорение грузика пропорциональны соответственно угловой скорости и ускорению (то есть $d θ / d t$ и $d^{2} θ / d t^{2}$ ), а также длине маятника. Действительно, если увеличить длину вдвое, при такой же угловой скорости реальная скорость грузика увеличится также вдвое. Это означает, что реальное ускорение равно $l \cdot d^{2} θ / d θ^{2}$ и второй закон Ньютона в проекции на направление, перпендикулярное к стержню, запишется в следующем виде:

m l \frac{d^{2} θ}{d t^{2}} = - m g sin θ

Знак минус в правой части связан с тем фактом, что сила направлена в направлении уменьшения

θ

. Сокращая на

m

и перенося

l

в правую часть, имеем:

\frac{d^{2} θ}{d t^{2}} = \frac{g}{l} sin θ

Будем считать для простоты, что у нас

l = g

\frac{d^{2} θ}{d t^{2}} = - sin θ

Потенциальная энергия имеет вид:

U (θ) = - \int sin θ d θ = - cos θ

Полная энергия:

E (θ, v) = \frac{v^{2}}{2} - cos θ .

Здесь мы считаем, что

m = 1

Построим график полной энергии и её линии уровня.

import plotly
import plotly.graph_objs as go
from plotly.offline import iplot as plot
from plotly.offline import init_notebook_mode

init_notebook_mode()

import numpy as np

x, y = np.mgrid[-2*np.pi:2*np.pi:300j, -2:2:300j]
surface = go.Surface(
    x=x, y=y, z=-np.cos(x)+y**2/2, 
    contours=dict(
        y=dict(show=False), 
        x=dict(show=True), 
        z=dict(show=True, project=dict(z=True)))
)
plot([surface])

Построим также график функции потенциальной энергии и линии уровня полной энергии один под другим (они связаны).

import matplotlib.pyplot as plt
import numpy as np
import qqmbr.odebook as ob
# see https://github.com/ischurov/qqmbr/blob/master/qqmbr/odebook.py

fig = plt.figure(figsize=(6,12))
ax=fig.add_subplot(211) #FIGURE 1
plt.plot([3.14,3.14],[4,-4],'--k',[-3.14,-3.14],[4,-4],'--k',[0,0],[4,-4],'--k')
ob.axes4x4(labels=('\\theta', 'U'))
ob.mplot(np.linspace(-4, 4), lambda x: -np.cos(x))
ax=fig.add_subplot(212) #FIGURE 2
ob.axes4x4(labels=('\\theta','v'))
ob.mcontour(np.linspace(-5, 5, 300), np.linspace(-5, 5, 300),
    lambda x, y: y**2/2-np.cos(x), 
    levels=list(np.arange(-1, 5, 0.5)) + [0.1, -1 + 0.1, -1 + 0.01])
plt.plot([3.14,3.14],[4,-4],'--k',[-3.14,-3.14],[4,-4],'--k',[0,0],[4,-4],'--k')

Рис. 7.7: Потенциальная энергия и линии уровня полной энергии

Вопрос 1. В какую сторону должны быть направлены стрелки на картинке в фазовом пространстве?

7.3.1Интерпретация

Положениями равновесия системы Ньютона являются точки в фазовом пространстве, в которых сила равна нулю и скорость равна нулю (то есть это точки горизонтальной оси, в которых

F (x) = 0

). Они же соответствуют точкам максимума и минимума для потенциальной энергии (что неудивительно, потому что последняя с точностью до знака равна интегралу от силы). Вблизи точек минимума потенциальной энергии система ведёт себя подобно осциллятору (см. пример 1), а вблизи точек максимума потенциальной энергии — подобно системе из примера 2.

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

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

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

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

В окрестности точки $(0, 0)$ (точнее, любой точки $(2 π k, 0)$ имеем минимум потенциальной энергии и точку устойчивого равновесия (фазовые кривые близки к окружностям), в окрестности $(π + 2 π k, 0)$ — максимум потенциальной энергии и точку неустойчивого равновесия (фазовые кривые близки к гиперболам). Фазовые кривые, похожие на окружности, соответствуют колебательным движениям маятника. Фазовые кривые, похожие на синусоиды, соответствуют вращению маятника («солнышко» на качелях).

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

Вы можете поэкспериментировать с математическим маятником с помощью интерактивной демонстрации.

← Предыдущая глава Следующая глава →

Обыкновенные дифференциальные уравнения Интерактивный учебник

7Консервативные системы с одной степенью свободы

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

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

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

7.3.1Интерпретация

Обыкновенные дифференциальные уравнения
Интерактивный учебник