Побудувати математичну модель, яка описує рух підвісного крану по жорсткозакріпленій підвісній рейці. Данна модель повинна дозволяти визначити положення підвісу та вантажу в довільний момент часу.
Вхідні параметри:
Задача – знайти ці невідомі функції $x_1(t)$, $x_2(t)$, $y_2(t)$, які описують коливання математичного маятника з рухомою точкою підвісу.
Де:
Для отримання рівняння руху для розгойдуємо системи використовуємо лагранжевий формалізм.
Складаємо Лагранжіан системи: $[Lg = T-U]$, де $T$ – кінетична енергія системи, а $U$ – потенціальна енергія системи.
\begin{equation} \begin{matrix} T = T_1+T_2\\ T_1 = \frac{m_1V_1^2}{2}=\frac{m_1\dot{x}^2_1}{2}\\ T_2 = \frac{m_2V_2^2}{2}=\frac{m_2(\dot{x}^2_2+\dot{y}^2_2)}{2}\\ U = U_1+U_2\\ U_1 = 0\\ U_2 = -m_2gy_2\\ \dot{x}_2=\dot{x}_1-l\dot{\alpha}cos(\alpha)\\ \dot{y}_2=-l\dot{\alpha}sin(\alpha)\\ Lg = \frac{m_1\dot{x}_1^2}{2}+\frac{m_2}{2}(\dot{x}_1^2+l^2\dot{\alpha}^2+2l\dot{x}\dot{\alpha}cos(\alpha))+m_2glcos(\alpha)\\ \end{matrix} \end{equation}Складемо рівняння руху:
\begin{equation} \begin{matrix} \frac{d}{dt}\left(\frac{\partial Lg}{\partial\dot{\alpha}}\right)-\frac{\partial Lg}{\partial\alpha}=0\\ \frac{d}{dt}\left(\frac{\partial Lg}{\partial\dot{x}_1}\right)-\frac{\partial Lg}{\partial x_1}=0\\ \frac{\partial Lg}{\partial\alpha}=m_2l^2\dot{\alpha}+m_2l\dot{x}cos(\alpha)\\ \frac{d}{dt}\left(\frac{\partial Lg}{\partial\dot{\alpha}}\right)=m_2l^2\ddot{\alpha}-m_2l\ddot{x}cos(\alpha)-m_2l\dot{x}\dot{\alpha}sin(\alpha)\\ \frac{\partial Lg}{\partial t}=-m_2l\dot{x}\dot{\alpha}sin(\alpha)-m_2glsin(\alpha)\\ m_2l^2\ddot{\alpha}-m_2l\ddot{x}cos(\alpha)+m_2glsin(\alpha)&(1)\\ \frac{\partial Lg}{\partial x_1}=m_1\dot{x}_1+m_2\dot{x}_1+m_2l\dot{\alpha}cos(\alpha)\\ \frac{d}{dt}\left(\frac{\partial Lg}{\partial\dot{x}_1}\right)=(m_1+m_2)\ddot{x}_1+m_2l\ddot{\alpha}cos(\alpha)-m_2l\ddot{\alpha}^2sin(\alpha)\\ (m_1+m_2)\ddot{x}_2+m_2l\ddot{\alpha}cos(\alpha)-m_2l\dot{\alpha}^2sin(\alpha)=0&(2)\\ \end{matrix} \end{equation}В результаті підстановки задача представляє собой систему двох диференціальних рівнянь 2-го порядку, які є нелінійними. Точного елементарного розв'язку данна система не має.
Визначимо початкові умови:
\begin{equation} \begin{matrix} & \left\{ \begin{matrix} x_1(0)=x_0\\ \alpha(0)=\alpha_0\\ \dot{x}_1(0)=v_0\\ \dot{\alpha}(0)=\omega_0 \end{matrix} \right. \end{matrix} (3) \end{equation}Постановка (1), (2) і (3) описує поведінку побудованної системи у довільний момент часу $t>0$.
Розв'яжемо отриману систему:
\begin{equation} \begin{matrix} & \left\{ \begin{matrix} \frac{d^2}{dt^2}x(t)+l\frac{d^2}{dt^2}\alpha(t)+g\alpha(t)=0\\ \frac{d^2}{dt^t}x(t)(m_1+m_2)+lm_2\frac{d^2}{dt^2}\alpha(t)=0\\ x_1(0)=x_0\\ \alpha(0)=\alpha_0\\ \dot{x}_1(0)=v_0\\ \dot{\alpha}(0)=\omega_0 \end{matrix} \right. \end{matrix} \end{equation}Отримаємо:
\begin{equation} \begin{matrix} &\alpha(t)=\frac{v_0\sqrt{m_1}\sqrt{l}sin\left(\frac{\sqrt{g}\sqrt{m_1+m_2}t}{\sqrt{m_1}\sqrt{l}}\right)}{\sqrt{g}\sqrt{m_1+m_2}}+a_0cos\left(\frac{\sqrt{g}\sqrt{m_1+m_2}t}{\sqrt{m_1}\sqrt{l}}\right)\\ &x(t)=-\frac{1}{m_1+m_2}\left(\frac{sin\left(\frac{\sqrt{g}\sqrt{m_1+m_2}t}{\sqrt{m_1}\sqrt{l}}\right)l^\frac{3}{2}v_0\sqrt{m_1}m_2}{\sqrt{g}\sqrt{m_1+m_2}}+cos\left(\frac{\sqrt{g}\sqrt{m_1+m_2}t}{\sqrt{m_1}\sqrt{l}}\right)la_0m_2-\frac{lv_0m_2tm_1}{m_1+m_2}-\frac{lv_0m^2_2t}{m_1+m_2}-\frac{(a_0lm_2+m_1x_0+m_2x_0)m_1}{m_1+m_2}-\frac{(a_0lm_2+m_1x_0+m_2x_0)m_2}{m_1+m_2}\right)\\ \end{matrix} \end{equation}Знайшовши розв'язок, оголошуємо основні функції a(t) та x(t):
In [1]:
# Імпортуємо необхідні модулі.
from ipywidgets import *
from numpy import sin, cos, sqrt, pi, radians, arange
import matplotlib.pyplot as plt
from matplotlib import rc
font = {'family': 'Verdana',
'weight': 'normal'}
rc('font', **font)
# Оголошуємо функції.
def ar():
return v0*sqrt(m1)*sqrt(l)*sin(sqrt(g)*sqrt(m1+m2)*t/(sqrt(m1)*sqrt(l)))/(sqrt(g)*sqrt(m1+m2))+a0*cos(sqrt(g)*sqrt(m1+m2)*t/(sqrt(m1)*sqrt(l)))
def xr():
return -(1.0/(m1+m2))*(sin(sqrt(g)*sqrt(m1+m2)*t/(sqrt(m1)*sqrt(l)))*l**(3.0/2.0)*v0*sqrt(m1)*m2/(sqrt(g)*sqrt(m1+m2))+cos(sqrt(g)*sqrt(m1+m2)*t/(sqrt(m1)*sqrt(l)))*a0*l*m2-(l*v0*m2*t*m1/(m1+m2))-(l*v0*(m2**2.0)*t/(m1+m2))-((a0*l*m2+m1*x0+m2*x0)*m1/(m1+m2))-((a0*l*m2+m1*x0+m2*x0)*m2/(m1+m2)))
А також функції для визначення положення підвісу та вантажу у довільний момент часу:
In [2]:
def x1():
return xr()
def x2():
return x1()+l*sin(ar())
def y2():
return -l*cos(ar())
Визначаємо початкові умови:
In [3]:
g = 9.8 # Час прискорення вільного падіння.
m1 = 3.0 # Маса підвісу.
m2 = 2.0 # Маса вантажу.
l = 6.0 # Довжина тросу підвісу.
a0 = pi/4.0 # Кут початкового відхилення від положення рівноваги.
v0 = 2.0 # Початкова швидкість.
x0 = 0.0 # Початкова координата по вісі Ox.
global_length, delta = 15.0, 0.01 # global_length - загальний час руху, delta - крок.
t = arange(x0, global_length, delta) # t - інтервал часу від 'x0' до 'global_length' з кроком 'delta'.
Будуємо графіки:
In [4]:
figure1 = plt.figure()
plot1 = figure1.add_subplot(3, 2, 1)
plot2 = figure1.add_subplot(3, 2, 2)
plot3 = figure1.add_subplot(3, 1, 2)
for ax in figure1.axes:
ax.grid(True)
plt.subplots_adjust(top=2.0, right=2.0, wspace=0.10, hspace=0.25)
plot1.plot(t, xr(), 'g')
plot2.plot(t, ar(), 'g')
plot3.plot(x1(), [0.01 for i in range(len(t))], 'b')
plot3.plot(x2(), y2(), 'r')
plot1.set_title(u'xr(t)')
plot2.set_title(u'ar(t)')
plot3.set_title(u'Траекторія руху маятника')
plt.show()
Отримана система містить рівняння другого порядку. Отже, для побудови наближенного розв'язку методом Ейлера необхідно звести ці рівняння до першого порядку, додавши нові змінні:
\begin{equation} \begin{matrix} y=\dot{x}\\ \dot{y}=\ddot{x}\\ \beta=\dot{\alpha}\\ \dot{\beta}=\ddot{\alpha}\\ \end{matrix} \Rightarrow \begin{matrix} l\dot{\beta}+\dot{y}+g\alpha=0\\ (m_1+m_2)\dot{y}+m_1l\dot{\beta}=0\\ \end{matrix} \end{equation}\begin{equation} \begin{matrix} \dot{x}=y&(4)\\ \dot{\alpha}=\beta&(5)\\ (m_1+m_2)\dot{y}+m_1l\dot{\beta}=0\\ \left(\frac{m_1+m_2}{m_1l}\right)\dot{y}+\left(\frac{m_1l\dot{\beta}}{m_1l}\right)=0\\ \dot{\beta}=-\frac{m_1+m_2}{m_1l}\dot{y}&(6)\\ l\dot{\beta}+\dot{y}+g\alpha=0\\ l\left(-\frac{m_1+m_2}{m_1l}\dot{y}\right)+\dot{y}+g\alpha=0\\ -\frac{m_1}{m_2}+\dot{y}=-g\alpha\\ \dot{y}=\frac{m_2}{m_1}g\alpha&(7)\\ \end{matrix} \end{equation}Об'єднуємо у систему рівняння $(4)$, $(5)$, $(6)$, $(7)$ та визначаємо початкові умови:
\begin{equation} \left\{ \begin{matrix} \dot{x}=y\\ \dot{\alpha}=\beta\\ \dot{y}=\frac{m_2}{m_1}g\alpha\\ \dot{\beta}=-\frac{m_1+m_2}{m_1l}\frac{g}{l}\alpha\\ \end{matrix} \right. \begin{matrix} x(0)=x_0\\ \alpha(0)=\alpha_0\\ y(0)=v_0\\ \beta(0)=u_0\\ \end{matrix} \end{equation}Задача - знайти розв'язок деякої функції $\dot{y}(t)=f(t,y(t))$, де $t(t_0)=y_0$, на інтервалі від $x_0$ до $n$. Для цього необхідно побудувати цикл, кожна ітерація якого, залежить від результату попередньої: $y_{i+1}=y_i+\Delta tf(t_i, y_i)$, де $\Delta t$ - крок ітерації.
\begin{equation} \left\{ \begin{matrix} x_{i+1}=x_i+\Delta ty_i\\ \alpha_{i+1}=\alpha_i+\Delta t\beta_i\\ y_{i+1}=y_i+\Delta t\frac{m_2}{m_1}g\alpha_i\\ \beta_{i+1}=\beta_i-\Delta t\frac{m_1+m_2}{m_1}\frac{g}{l}\alpha_i\\ \end{matrix} \right. \end{equation}
In [5]:
x = arange(x0, global_length, delta)
a = arange(x0, global_length, delta)
y = arange(x0, global_length, delta)
b = arange(x0, global_length, delta)
x[0], a[0], y[0], b[0] = x0, a0, v0, sqrt(g/l)
for i in range(0, len(t)-1):
x[i+1] = x[i] + delta * y[i]
a[i+1] = a[i] + delta * b[i]
y[i+1] = y[i] + delta * (m2 / m1) * g * a[i]
b[i+1] = b[i] - delta * ((m1 + m2) / m1) * (g/l) * a[i]
def aproxy_x1():
return x
def aproxy_x2():
return x+l*sin(a)
def aproxy_y2():
return -l*cos(a)
Будуємо графіки:
In [6]:
figure2 = plt.figure()
plot4 = figure2.add_subplot(3, 2, 1)
plot5 = figure2.add_subplot(3, 2, 2)
plot6 = figure2.add_subplot(3, 1, 2)
for ax in figure2.axes:
ax.grid(True)
plt.subplots_adjust(top=2.0, right=2.0, wspace=0.10, hspace=0.25)
plot4.plot(t, xr(), 'r')
plot4.plot(t, x, 'b')
plot5.plot(t, ar(), 'r')
plot5.plot(t, a, 'b')
plot6.plot(x1(), [0.01 for i in range(len(t))], 'r')
plot6.plot(aproxy_x1(), [0.01 for i in range(len(t))], 'b')
plot6.plot(x2(), y2(), 'r', label=u'Точний розв\'язок')
plot6.plot(aproxy_x2(), aproxy_y2(), 'b', label=u'Метод Ейлера')
plot4.set_title(u'xr(t)')
plot5.set_title(u'ar(t)')
plot6.set_title(u'Порівняння точноі та наближеної траекторії руху маятника')
plot6.legend(loc='upper right')
plt.show()