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])

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

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

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

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

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

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