Математичне моделювання маятника


Змістовна постановка

Побудувати математичну модель, яка описує рух підвісного крану по жорсткозакріпленій підвісній рейці. Данна модель повинна дозволяти визначити положення підвісу та вантажу в довільний момент часу.

Вхідні параметри:

  • $l$ – довжина тросу підвісу
  • $m_1$ – масса підвісу
  • $m_2$ – масса вантажу

Концептуальна постановка

  1. Трос будемо вважати невагомим неростяжним стержнем, жорсткозакріпленим у точці ванатжу та у точці підвісу.
  2. Нехтуємо силами опору тертя у точці підвісу та силами опору зовнішнього середовища.
  3. Модель будемо будувати у рамках класичної аналітичної механіки.

Математична постановка

Задача – знайти ці невідомі функції $x_1(t)$, $x_2(t)$, $y_2(t)$, які описують коливання математичного маятника з рухомою точкою підвісу.

Де:

  • $\alpha(t)$ – кут відхилення від положення рівноваги
  • $x_2(t) = x_1(t)+lsin(\alpha(t))$
  • $y_2(t) = -lcos(\alpha(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)))

А також функції для визначення положення підвісу та вантажу у довільний момент часу:

  • $x_1(t) = x(t)$
  • $x_2(t) = x_1(t)+lsin(\alpha(t))$
  • $y_2(t) = -lcos(\alpha(t))$

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} \left\{ \begin{matrix} l\ddot{\alpha}+\ddot{x}+g\alpha=0\\ (m_1+m_2)\ddot{x}+m_1l\ddot{\alpha}=0\\ \end{matrix} \right. \end{equation}

Отримана система містить рівняння другого порядку. Отже, для побудови наближенного розв'язку методом Ейлера необхідно звести ці рівняння до першого порядку, додавши нові змінні:

\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()