Прикладные дифференциальные уравнения

Домашнее задание №7 (бонусное)

Илья Щуров

Факультет компьютерных наук, Прикладная математики и информатика, 2021-22 учебный год

Страница курса

Задача 1: двойной маятник

Сделайте симуляцию двойного маятника. Положение маятника в фиксированный момент времени задаётся двумя числами: $\theta_1$ и $\theta_2$ — это углы наклона первой и второй компоненты маятника относительно вертикали. Вместо скоростей мы будем использовать так называемые обобщённые импульсы $p_{\theta_1}$ и $p_{\theta_2}$. Что это такое — сейчас не очень важно; мы не будем выводить следующие формулы движения, а просто поверим Википедии. (Скажу только, что есть такой прикольный способ отыскивать уравнение движения механических систем — через лагранжиан — но его обсуждение выходит далеко за рамки нашего курса.) Итак, формулы движения выглядят так:

$$ \begin{cases} {\dot \theta_1} = \dfrac{6}{m\ell^2} \dfrac{ 2 p_{\theta_1} - 3 \cos(\theta_1-\theta_2) p_{\theta_2}}{16 - 9 \cos^2(\theta_1-\theta_2)}\\ \\ {\dot \theta_2} = \dfrac{6}{m\ell^2} \dfrac{ 8 p_{\theta_2} - 3 \cos(\theta_1-\theta_2) p_{\theta_1}}{16 - 9 \cos^2(\theta_1-\theta_2)}\\ {\dot p_{\theta_1}} = -\frac{1}{2} m \ell^2 \left [ {\dot \theta_1} {\dot \theta_2} \sin (\theta_1-\theta_2) + 3 \frac{g}{\ell} \sin \theta_1 \right ]\\ {\dot p_{\theta_2}} = -\frac{1}{2} m \ell^2 \left [ -{\dot \theta_1} {\dot \theta_2} \sin (\theta_1-\theta_2) + \frac{g}{\ell} \sin \theta_2 \right ]. \end{cases} $$

Здесь $m$ — масса каждого стержня, $\ell$ — длина каждого стержня, $g$ — ускорение свободного падения. Вы можете взять $m=1$ и $\ell = 1$ для определённости.

Обратите внимание, что в правой части уравнений на $\dot{p}_{\theta_1}$ и $\dot{p}_{\theta_2}$ встречаются $\dot \theta_1$ и $\dot \theta_2$. Их нужно просто подставить согласно первым двум уравнениям.

Часть 1

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

In [ ]:
from math import sin, cos
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from IPython.display import display, clear_output

%matplotlib inline
In [ ]:
l = 1
g = 9.8


def rhs(t, X):
    theta, theta_dot = X
    return np.array([theta_dot, -g / l * sin(theta)])
In [ ]:
def plot_pendulum(X, ax, color="C0"):
    theta, theta_dot = X
    ax.plot([0, l * sin(theta)], [0, -l * cos(theta)], "o-", color=color)
In [ ]:
# эта ячейка показывает анимацию в реальном времени
# она будет выполняться бесконечно долго
# остановите её с помощью Kernel → Interrupt

# вы можете использовать эту ячейку, чтобы сразу увидеть результат
# и иметь возможность быстро отлаживать вашу программу,
# а когда всё будет готово, сделать нормальную анимацию с помощью
# celluloid или matplotlib.animation, как вы это делали раньше

initial_condition = np.array([np.pi / 2, 0])
delta_t = 0.02
t = 0

fig, ax = plt.subplots(figsize=(10, 10))

system = ode(rhs)
system.set_initial_value(initial_condition)

while True:
    t += delta_t
    X = system.integrate(t)
    ax.cla()
    plot_pendulum(X, ax)
    ax.set_aspect("equal", "box")
    ax.set(xlim=(-3, 3), ylim=(-3, 3))
    display(fig)
    clear_output(wait=True)
    # plt.pause(0.02)

Часть 2

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

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

Задача 2: маятник с трением

Рассмотрим математический маятник с трением. Его динамика задаётся уравнением $$\ddot \theta = -\sin \theta - \mu \dot \theta,$$ где $\mu > 0$ — коэффициент трения. (Мы опять для простоты рассматриваем вязкое трение, сила которого пропорциональна скорости.)

Часть 1

Возьмите $\mu=0.2$. Нарисуйте с помощью компьютера фазовый портрет системы в области $\theta \in (-8, 8)$, $\dot \theta \in (-5, 5)$. Нарисуйте на этом фазовом портрете траекторию с начальным условием $\theta_0=-7.5$, $\dot \theta_0 = 4$. Опишите словами, что происходит с маятником, когда он движется по этой траектории.

Часть 2

Полная механическая энергия системы равна $$E(\theta, \dot \theta) = \frac{\dot \theta^2}{2} - \cos \theta.$$ Докажите аналитически, что если система не находится в положении равновесия, то с течением времени полная механическая энергия строго убывает при любом $\mu > 0$. (Подсказка: вам понадобится производная функции вдоль векторного поля.) Из этого факта можно вывести, что все решения стремятся к положениям равновесия (но этого сейчас делать не нужно).

Часть 3

Найдите все особые точки (положения равновесия) системы, лежащие в области $\theta \in [0, 2\pi)$. Каким положениям маятника соответствуют эти особые точки? Найдите линеаризацию системы в каждой особой точке. Определите тип каждой особой точки в зависимости от $\mu > 0$. При каком значении $\mu$ маятник совершает бесконечно много колебаний при стремлении к устойчивому равновесию?

Часть 4

Постройте фазовые портреты линеаризованных систем в окрестностях особых точек, найденных в прошлой части, на той же картинке, на которой построен фазовый портрет исходной системы для $\mu=0.2$. (Особая точка линеаризованной системы должна совпадать с особой точкой исходной — сдвиньте линеаризованную систему путём подходящей замены координат если это необходимо.) Для особой точки, являющейся седлом, нарисуйте собственные векторы линеаризации (отложите их от этой особой точки).

Часть 5

Найдите и нарисуйте с помощью компьютера (например, с помощью бисекции отрезка) сепаратрисы седла для исходной системы (c $\mu=0.2$). Убедитесь (по картинке), что сепаратрисы касаются соответствующих собственных векторов.

Часть 6

Верно ли, что если в начальный момент времени маятник не находится в положении равновесия, то со временем он обязательно будет стремиться к положению равновесия, при котором он висит вертикально вниз? Ответ обосновать, опираясь на результаты предыдущих частей.

In [ ]: