«Физика для программистов» — как физтехи применяют её в приложениях. Маятники

Простой

8 мин

1.6K

Аннотация

Данная статья входит в цикл, освещающий задачи на моделирование физических процессов на факультете МТФИ ВШПИ. В этой части речь пойдёт про задачу моделирования поведения маятника: коротко разберём теорию, которая лежит в основе модели, немного подумаем над архитектурой и напишем небольшое приложение на связке Python + Tkinter. Реализация будет поддерживать исследование различных маятников с помощью самописных динамических графиков, в которые пользователь может ввести собственные формулы.

Введение

Эта формулировка была предложена в качестве одной из задач для моделирования на первом семестре курса «Цифровизация физических процессов» на факультете ВШПИ МФТИ:

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

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

Теория

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

Для дальнейших расчётов введём обозначения:
– момент инерции маятника относительно оси,
– расстояние от центра масс до оси,
– масса маятника,
$\phi$ – угол отклонения маятника от вертикальной оси,
$\ddot\phi$ – вторая производная угла по времени (угловое ускорение),
– ускорение свободного падения.

Тогда из уравнения динамики твёрдого тела получаем:

$I\ddot\phi = -mglsin\phi \rightarrow \ddot\phi = -\frac{mgl}{I}sin\phi = -\omega^{2}sin\phi$ (1)

Принято обозначать $\omega^{2} = \frac{mgl}{I}$ .

Эта почти та формула, которую мы будем использовать для моделирования, чуть позже мы её дополним, однако в теории чаще используются другие формулы:
$\phi=\Phi_{0}cos(\omega t + \phi_{0})$ (2)
$\phi=\Phi_{0}sin(\omega t + \phi_{0}')$
$\phi=Acos(\omega t) + Bsin(\omega t)$ ,
где, грубо говоря, $\Phi_{0}$ отвечает за максимальное отклонение, а $\phi_{0}$ за начальную позицию внутри цикла.
Легко видеть, что все эти формулы — решения описанного выше дифференциального уравнения при $\phi$ ≪ 1 радиана $\rightarrow sin\phi \approx \phi$ . Действительно, такое допущение, позволяет элегантно разрешить полученное дифференциальное уравнение аналитически, получив описанные выше формулы. Можете проверить их, продифференцировав дважды.

Наиболее нетривиальный параметр — это момент инерции. Напомню, что

$I = \int_{m}{r^{2}dm}$ .

Например, для математического маятника

$I = mr^{2}$ , т.к. вся масса удалена на от оси. Тогда по (2):

$\phi=\Phi_{0}cos(\sqrt{\frac{mgl}{mr^{2}}}t) = |r = l| = \Phi_{0}cos(\sqrt{\frac{g}{l}}t)$ (3),

отсюда легко получить период колебаний:

$T = 2\pi\sqrt{\frac{l}{g}}$ ,

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

В исходной задаче осталось учесть только влияние среды. В модели вязкого трения получаем:

$\ddot\phi = -\omega^{2}sin\phi - \gamma\dot\phi$ (4),

где $\gamma$ - коэффициент вязкого трения. Положив $\gamma$ = 0 мы получим (1), так что будем моделировать сразу общий случай по формуле (4).

Моделирование

Выбор технологий

Желания возится с точностью не было, как и необходимости в чрезмерной производительности, так что в качестве базы был выбран Python. За визуализацию будет отвечать Tkinter, здесь строго субъективно, я им никогда не пользовался, так что хотелось попробовать.

Архитектура

Так как мы пишем простенькое приложение, которое всё держит в рантайме, то ничего сильно выдумывать не стоит. Очевидно, что потоки отрисовки (view) и модели (model) стоит разделить. Действительно, представим, что поток отрисовки буксует, тогда модель будет замедляться, а хотелось бы, чтобы модель успевала обрабатывать необходимое число шагов в реальном времени, сохраняя при этом выбранный пользователем уровень дискретизации (в адекватных пределах). Пусть мы разделили потоки, на каждом шаге view будет рисовать актуальные данные, даже если она пропустила некоторые шаги обновления model, не замедляя его. Кроме того, при такой архитектуре мест, где могут проявиться проблемы с concurrency фактически нет. Для создания нескольких потоков использовалась библиотека threading.

Итого, наше приложение — это 2 потока и, соответственно, 2 движка, view и model, каждый шаг оба движка стараются обновить всех зарегистрированных entity. Фактически мы будем реализовывать MVC, с учётом того, что View-Controller представляет Tkinter.

Модель

Если отбросить довольно душную отрисовку и тривиальную абстракцию движков, то каждый маятник выполняет каждый тик одну из этих функций в зависимости от его типа:

def step_real(self):
    o2 = self.main.params["g"] * self.main.mass * self.main.mass_center_remoteness * self.main.size / self.main.inertia_moment
    self.main.angle_acceleration = - o2 * sin(self.main.angle) - self.main.gamma * self.main
    self.main.angle_velocity += self.main.angle_acceleration * self.main.params["dt"]
    self.main.angle += self.main.angle_velocity * self.main.params["dt"]
    
def step_theory(self):
    o2 = self.main.params["g"] * self.main.mass * self.main.mass_center_remoteness * self.main.size / self.main.inertia_moment
    self.main.angle = self.main.start_angle * cos(sqrt(o2) * (self.main.params["model_tick"] - self.main.start_tick) * self.main.params["dt"])

Это те самые 2 формулы (4) и (2) из теоретических знаний, которые были разобраны выше.

UI

Задача

Итак, у нас есть две формулы, мы хотим уметь создавать/удалять/настраивать маятники, менять параметры окружения, проводить симуляцию, получать и анализировать данные. Разделим наше приложение на 4 части: workspace по центру и 3 контрольные панели по разным сторонам: снизу, справа и слева. Справа будут настройки выбранного маятника, слева — работа с данными, а снизу — параметры окружения.

Workspace

На workspace располагаются маятники. Workspace обладает функционалом создавать, перетаскивать маятники и минимальной возможностью настройки, также на workspace должны отображаться названия маятников, чтобы их можно было отличать. Маятник представляет из себя 6 частей: ось, сам маятник, противовес, центр масс, конец маятника и имя. При этом положения маятника, противовеса, центра масс и конца маятника необходимо менять во время отрисовки на актуальные данные от model.

Настройки маятника

Все основные параметры, которые есть у маятника, можно менять, в них входит: длина, расстояние до оси, масса и начальный угол. Кроме того, по умолчанию стоит «полная симуляция», но если хочется сравнить с теоретической моделью, этот флажок можно выключить. В этом же окне располагаются основные значения состояния каждого из маятников. Важно отметить, что рядом с каждым элементом указан код, по которому их можно получить, это будет важно в дальнейшем для написания собственных формул. В самом низу можно применить изменения, перезапустить маятник или удалить его с workspace.

Управление окружением

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

Поля формульного ввода

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

Итак, формула должна иметь вид выполняемого в Python математического выражения, где в качестве операндов могут выступать записи в виде @{id}.{поле}, ну и, конечно же, числа. Доступны следующие операции: «+», «-», «·», «/», «^», также доступны скобки. В целом, с небольшой модификацией кода можно добавить любые другие функции, которые поддерживает Python. Коротко говоря, программа парсит такие выражения и вызывает функцию eval() для расчёта значения. Если ввести некорректную формулу, поле ввода станет красным, сигнализируя об ошибке.

Подсчёт значения параметра вида @{id}.{field} происходит так:

def calculate_value(self, value_code):
    if is_number(value_code):
        return True, float(value_code)
    try:
        entity, target_name = value_code.split(".")
    except ValueError:
        return False, 0
    if not (entity in self.params["model"].entities and target_name in self.params["model"].entities[entity].targets):
        return False, 0
    return True, self.params["model"].entities[entity].targets[target_name].value

Соответственно, полный парсер выглядит так: (tokens = ["+", "-", "*", "/", " ", "(", ")"])

def calculate_target(self, target):
    if target == "":
        return False, 0
    value = ""
    command = ""
    for ch in target:
        is_token = False
        for token in self.tokens:
            if ch == token:
                if value != "":
                    ok, result = self.calculate_value(value)
                    if not ok:
                        return False, 0
                    command += "(" + str(result) + ")"
                value = ""
                command += str(ch)
                is_token = True
                break
        if not is_token:
            value += ch
    if value != "":
        ok, result = self.calculate_value(value)
        if not ok:
            return False, 0
        if result > 0:
            command += str(result)
        else:
            command += "(" + str(result) + ")"
    try:
        target_value = eval(command)
        return True, target_value
    except:
          return False, 0

Мы идём по каждому символу и ищем записи вида @{id}.{field}.

Графики

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

Классическая часть

В основном в статьях разбирается пример построения статических графиков, которые в целом можно просто вставить на canvas из matplotlib, однако основные идеи построения графика, такие как построение осей, делений и т. п., вы можете прочитать в них, я же хотел бы заострить внимание на некоторых интересных фичах, которые я применил в своей реализации.

Автофокус на графике

Чтобы удерживать график в удобном для пользователя формате, достаточно двух параметров масштаба и линии центра. Проанализировав особенности данных, которые теоретически можно получить, было принято решение сделать 4 основных режима отображения: удерживать 0 в центре, удерживать 0 сверху/снизу и не удерживать 0.

Ниже представлен код обновления y_scale и draw_axis_y — полей, отвечающих за масштаб графика и линию центра для разных режимов отображения (для более общей настройки y_scale считается в экспоненциальном виде):

def center_behaviour(self, max_y, min_y):
    self.draw_axis_y = 0
    if max_y != 0 and min_y != 0:
        self.y_scale = log2(min(abs(self.canvas_size_y / 2 / max_y), abs(self.canvas_size_y / 2 / min_y)) * 2 / 3)

def axis_up_behaviour(self, max_y, min_y):
    self.draw_axis_y = max_y / 2
    if max_y != 0:
      self.y_scale = log2(abs(self.canvas_size_y / max_y) * 2 / 3)

def axis_down_behaviour(self, max_y, min_y):
    self.draw_axis_y = min_y / 2
    if min_y != 0:
      self.y_scale = log2(abs(self.canvas_size_y / min_y) * 2 / 3)

def custom_behaviour(self, max_y, min_y):
    self.draw_axis_y = (min_y + max_y) / 2
    if min_y != max_y:
        self.y_scale = log2(abs(self.canvas_size_y / abs(min_y - max_y)) * 2 / 3)

Каждый из режимов удобен в своих условиях, а в совокупности они успешно покрывают основные потребности.

Получившийся график

Пример графика, который выводит значения угла отклонения маятника с id=p2

Чтобы избежать информационной перегруженности графика, из делений было решено оставить только ноль и минимум/максимум на canvas'е, по этим данным уже можно что-то быстро прикинуть. Точные данные в каждой точки можно посмотреть, если навестись на нужное место на canvas'е, как и в matplotlib. В целом, компонент графика получился функциональным, наглядным и удобным для управления. Конечно, в нём есть что оптимизировать и доделывать, но в контексте приложения этого достаточно.

Заключение

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

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

Ссылки

GitHub проекта

Теги:

Хабы: