Введение
Двойной маятник - это классическая задача в физике и математике. Она включает два маятника, соединенных друг с другом, и движение второго маятника зависит от движения первого маятника. В этом лабораторном занятии мы будем использовать Python и Matplotlib для моделирования движения двойного маятника.
Советы по работе с ВМ
После запуска ВМ нажмите в верхнем левом углу, чтобы переключиться на вкладку Notebook и получить доступ к Jupyter Notebook для практики.
Иногда вам может потребоваться подождать несколько секунд, пока Jupyter Notebook не загрузится полностью. Валидация операций не может быть автоматизирована из-за ограничений Jupyter Notebook.
Если вы сталкиваетесь с проблемами во время обучения, не стесняйтесь обращаться к Labby. Оставьте отзыв после занятия, и мы оперативно решим проблему для вас.
Импортируем необходимые библиотеки
Начнем с импорта необходимых библиотек для нашей моделирования.
from collections import deque
import matplotlib.pyplot as plt
import numpy as np
from numpy import cos, sin
import matplotlib.animation as animation
Задаем параметры
Далее мы определим параметры для нашей моделирования. Эти параметры включают ускорение свободного падения, длину и массу каждого маятника, а также интервал времени для моделирования.
G = 9.8 ## ускорение свободного падения, м/с^2
L1 = 1.0 ## длина первого маятника, м
L2 = 1.0 ## длина второго маятника, м
L = L1 + L2 ## максимальная длина объединенного маятника
M1 = 1.0 ## масса первого маятника, кг
M2 = 1.0 ## масса второго маятника, кг
t_stop = 2.5 ## на сколько секунд моделировать
history_len = 500 ## количество точек траектории для отображения
Определяем функцию для вычисления производных
Мы определим функцию, которая вычисляет производные системы в любое заданное время. Эта функция принимает текущее состояние системы (углы и угловые скорости каждого маятника) и возвращает производные этих значений.
def derivs(t, state):
dydx = np.zeros_like(state)
dydx[0] = state[1]
delta = state[2] - state[0]
den1 = (M1+M2) * L1 - M2 * L1 * cos(delta) * cos(delta)
dydx[1] = ((M2 * L1 * state[1] * state[1] * sin(delta) * cos(delta)
+ M2 * G * sin(state[2]) * cos(delta)
+ M2 * L2 * state[3] * state[3] * sin(delta)
- (M1+M2) * G * sin(state[0]))
/ den1)
dydx[2] = state[3]
den2 = (L2/L1) * den1
dydx[3] = ((- M2 * L2 * state[3] * state[3] * sin(delta) * cos(delta)
+ (M1+M2) * G * sin(state[0]) * cos(delta)
- (M1+M2) * L1 * state[1] * state[1] * sin(delta)
- (M1+M2) * G * sin(state[2]))
/ den2)
return dydx
Задаем начальные условия и интегрируем дифференциальное уравнение
Теперь мы зададим начальные условия для нашей моделирования. Мы зададим начальные углы и угловые скорости каждого маятника, а также интервал времени для моделирования. Затем мы интегрируем дифференциальное уравнение методом Эйлера, чтобы получить положение и скорость каждого маятника на каждом временном шаге.
## create a time array from 0..t_stop sampled at 0.02 second steps
dt = 0.01
t = np.arange(0, t_stop, dt)
## th1 and th2 are the initial angles (degrees)
## w10 and w20 are the initial angular velocities (degrees per second)
th1 = 120.0
w1 = 0.0
th2 = -10.0
w2 = 0.0
## initial state
state = np.radians([th1, w1, th2, w2])
## integrate the ODE using Euler's method
y = np.empty((len(t), 4))
y[0] = state
for i in range(1, len(t)):
y[i] = y[i - 1] + derivs(t[i - 1], y[i - 1]) * dt
Вычисляем положение каждого маятника
Теперь мы будем использовать положение и скорость каждого маятника на каждом временном шаге для вычисления координат x и y каждого маятника.
x1 = L1*sin(y[:, 0])
y1 = -L1*cos(y[:, 0])
x2 = L2*sin(y[:, 2]) + x1
y2 = -L2*cos(y[:, 2]) + y1
Настраиваем график
Теперь мы настроим график для нашей моделирования. Мы создадим фигуру с ограничениями по оси x и y, равными максимальной длине маятника, установим соотношение сторон равным и добавим сетку.
fig = plt.figure(figsize=(5, 4))
ax = fig.add_subplot(autoscale_on=False, xlim=(-L, L), ylim=(-L, 1.))
ax.set_aspect('equal')
ax.grid()
Определяем функцию для анимации
Мы определим функцию, которая будет анимировать движение двойного маятника. Эта функция будет принимать текущий временной шаг и использовать его для обновления положения каждого маятника. Она также обновит траекторию второго маятника.
line, = ax.plot([], [], 'o-', lw=2)
trace, = ax.plot([], [], '.-', lw=1, ms=2)
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
history_x, history_y = deque(maxlen=history_len), deque(maxlen=history_len)
def animate(i):
thisx = [0, x1[i], x2[i]]
thisy = [0, y1[i], y2[i]]
if i == 0:
history_x.clear()
history_y.clear()
history_x.appendleft(thisx[2])
history_y.appendleft(thisy[2])
line.set_data(thisx, thisy)
trace.set_data(history_x, history_y)
time_text.set_text(time_template % (i*dt))
return line, trace, time_text
Создаем анимацию
Теперь мы создадим анимацию с использованием функции FuncAnimation из Matplotlib.
ani = animation.FuncAnimation(
fig, animate, len(y), interval=dt*1000, blit=True)
plt.show()
Резюме
В этом лабораторном занятии мы использовали Python и Matplotlib для моделирования движения двойного маятника. Мы установили параметры для нашей моделирования, определили функцию для вычисления производных системы, интегрировали ОДУ методом Эйлера, вычислили положение каждого маятника на каждом временном шаге и создали анимацию движения двойного маятника.