Вывод и интегрирование уравнений движения эллиптического мятника в Python (использование sympy, scipy)

12 May 2024 in Python / Механика

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

Рассматриваемая механическая система имеет две степени свободы: её конфигурация определяется координатой $x$ тела 1 (брусок, движущийся поступательно вдоль оси $$Ox_0$) и углом отклонения подвеса тела 2 от вертикали $\phi$. Движение системы происходит в однородном поле силы тяжести.

Вывод уравнений движения

Подключение необходимых библиотек:

# Массивы и матрицы
import numpy as np
# Графики
import matplotlib.pyplot as plt
# Символьные вычисления
import sympy as sp
# Численное интегрирование дифференциальных уравнений
from scipy.integrate import solve_ivp

Параметры системы и обобщенные координаты

Уравнения движения системы запишем используя уравнения Лагранжа II-го рода. Для этого необходимо записать выражение для кинетической энергии системы $T(q_1,\dot q_1,q_2,\dot q_2)$, аналитически продифференцировать его по обобщенной скорости, времени, обобщенной координате: \[\frac{d}{dt} \frac{\partial T}{\partial \dot q_i} -\frac{\partial T}{\partial q_i} = Q_{i}, \quad i=1,2\]

Уравнения движения рассматриваемой механической системы будут выводиться с использованием библиотеки sympy. Все переменные и параметры, которые будут входить в уравнения движения, необходимо объявить до их использования при помощи функции var или symbols. Рассматриваемая система описывается следующими параметрами:

масса тела 1 ($m_1$);
масса тела 2 ($m_2$);
длина подвеса тела 2 ($l = AB$);
ускорение свободного падения ($g$).

Таким образом, объявление этих параметров-символов будет иметь вид:

sp.var('m1,m2,l,g')

Также необходимо объявить обобщенные координаты, которые являются функциями времени. Это делается при помощи функции Function:

x   = sp.Function('x')
phi = sp.Function('phi')

После объявления параметров и обобщенных координат (неизвестных пока функций времени) можно начать формирование выражений, необходимых для построения уравнений движения.

Кинематика (скорости тел)

Координатный столбец положения тела 1 в системе координат $Ox_0y_0$ в sympy определим при помощи объекта типа Матрица:

r1 = sp.Matrix([x(t),0])

При формировании этого выражения указывается, что обобщенная координата – объявленная ранее функция x является функцией времени.

Координатный столбец положения тела 2 в системе координат $Ox_0y_0$:

r2 = r1 + l*sp.Matrix([sp.sin(phi(t)),sp.cos(phi(t))])

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

v1 = sp.diff(r1,t) 
v1

\[\left[\begin{matrix}\frac{d}{d t} x{\left(t \right)}\\0\end{matrix}\right]\]

Аналогично определим координатный столбец вектора скорости тела 2:

v2 = sp.diff(r2,t)
v2

\[\displaystyle \left[\begin{matrix}l \cos{\left(\phi{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} + \frac{d}{d t} x{\left(t \right)}\\- l \sin{\left(\phi{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)}\end{matrix}\right]\]

Левая часть уравнений Лагранжа

Кинетическая энергия рассматриваемой системы определяется суммой кинетических энергий тела 1 и тела 2:

T = m1*v1.dot(v1)/2 + m2*v2.dot(v2)/2
T

\[\frac{m_{1} \left(\frac{d}{d t} x{\left(t \right)}\right)^{2}}{2} + \frac{m_{2} \left(l^{2} \sin^{2}{\left(\phi{\left(t \right)} \right)} \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2} + \left(l \cos{\left(\phi{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} + \frac{d}{d t} x{\left(t \right)}\right)^{2}\right)}{2}\]

Используя функцию diff найдем левую часть уравнений Лагранжа II-го рода для первой обобщенной координаты $\phi$: \[\frac{d}{dt} \frac{\partial T}{\partial \dot \phi} -\frac{\partial T}{\partial \phi} = Q_\phi\]

eq1 = sp.diff(sp.diff( T, sp.diff(phi(t),t) ),t) - sp.diff(T,phi(t))
# Упрощаем
eq1 = sp.simplify(eq1)
eq1

\[l m_{2} \left(l \frac{d^{2}}{d t^{2}} \phi{\left(t \right)} + \cos{\left(\phi{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} x{\left(t \right)}\right)\]

Левая часть уравнения Лагранжа II-го рода для второй обобщенной координаты $x$: \[\frac{d}{dt} \frac{\partial T}{\partial \dot x} -\frac{\partial T}{\partial x} = Q_x\]

eq2 = sp.simplify( sp.diff(sp.diff( T, sp.diff(x(t),t) ),t) - sp.diff(T,x(t)) )
eq2

\[l m_{2} \left(l \frac{d^{2}}{d t^{2}} \phi{\left(t \right)} + \cos{\left(\phi{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} x{\left(t \right)}\right)\]

Обобщенные силы

Механическая система в однородном поле силы тяжести, поэтому для определения обобщенных сил удобней использовать выражение для потенциальной энергии системы, определенной из условия того, что уровень y = 0 является уровнем нулевой потенциальной энергии: \[\Pi = - m_2 g l \cos \phi\]

P = -m2*g*r2[1]

Обобщенная сила для первой обобщенной координаты: \[Q_\phi = - \frac{\partial \Pi}{\partial \phi}\]

Q1 = -sp.diff(P,phi(t))

Обобщенная сила для второй обобщенной координаты: \[Q_x = - \frac{\partial \Pi}{\partial x}\]

Q2 = -sp.diff(P,x(t))

Приведение уравнений к форме Коши

Разрешим полученную систему дифференциальных уравнений (eq1 = Q1, eq2 = Q2) относительно старших производных, используя функцию solve:

сauchy_form = sp.simplify(sp.solve( (eq1 - Q1, eq2 - Q2), (sp.diff(phi(t),t,2), sp.diff(x(t),t,2) )))
сauchy_form

Функция solve возвращает решение в виде словаря: \[\left\{ \frac{d^{2}}{d t^{2}} \phi{\left(t \right)} : - \frac{g m_{1} \sin{\left(\phi{\left(t \right)} \right)}}{l m_{1} - l m_{2} \cos^{2}{\left(\phi{\left(t \right)} \right)} + l m_{2}} - \frac{g m_{2} \sin{\left(\phi{\left(t \right)} \right)}}{l m_{1} - l m_{2} \cos^{2}{\left(\phi{\left(t \right)} \right)} + l m_{2}} - \frac{l m_{2} \sin{\left(\phi{\left(t \right)} \right)} \cos{\left(\phi{\left(t \right)} \right)} \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2}}{l m_{1} - l m_{2} \cos^{2}{\left(\phi{\left(t \right)} \right)} + l m_{2}}, \ \frac{d^{2}}{d t^{2}} x{\left(t \right)} : \frac{g m_{2} \sin{\left(\phi{\left(t \right)} \right)} \cos{\left(\phi{\left(t \right)} \right)}}{m_{1} - m_{2} \cos^{2}{\left(\phi{\left(t \right)} \right)} + m_{2}} + \frac{l m_{2} \sin{\left(\phi{\left(t \right)} \right)} \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2}}{m_{1} - m_{2} \cos^{2}{\left(\phi{\left(t \right)} \right)} + m_{2}}\right\}\]

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

Сформируем список аргументов этой функции: первые 4 элемента этого списка – параметры системы, оставшиеся 4 элемента — обобщенные координаты и скорости (всего 8 параметров):

args = ( m1,m2,l,g,   phi(t),x(t),sp.diff(phi(t),t),sp.diff(x(t),t) )

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

accelerations = sp.lambdify( args, [ сauchy_form[sp.diff(phi(t),t,2)], сauchy_form[sp.diff(x(t),t,2)] ] )

Проверяем, как работает эта функция

accelerations(1.0, 3.0, 2, 9.81, 1.0,0, 0,0)

> [-5.284409988846435, 4.282768353189539]

Численное интегрирование

Объявим функцию правых частей для численного интегрирования системы двух дифференциальных уравнений второго порядка:

def dydt(t,q,p):            
    (d2phi, d2x) = accelerations(*(*p,*q))    
    return (*q[2:4], d2phi, d2x)  

Зададим параметры системы, начальные условия и проинтегрируем дифференциальные уравнения её движения:

# Параметры системы: m1, m2, l, g
p  = [3,1,2,9.81];
# Начальные условия phi0, x0, dphi/dt, dx/dt
q0 = [1.0, 0.0, 0, 0];

solution = solve_ivp(lambda t, q: dydt(t, q, p), [0, 10], q0, 
                     rtol = 1e-6, method="LSODA", t_eval = np.linspace(0,10,100))

Графики

plt.figure(figsize=[11,4])
plt.subplot(1,2,1)
plt.plot(solution.t,solution.y[0,:]*180/np.pi); plt.xlabel('t, c'); plt.ylabel('$\phi$, градус');
plt.subplot(1,2,2)
plt.plot(solution.t,solution.y[1,:]); plt.xlabel('t, c'); plt.ylabel('x, м');

Анимация

Дополнительные библиотеки для построения графических примитивов и анимации.

from matplotlib.animation import FuncAnimation
import matplotlib.lines as mlines
from matplotlib.patches import Wedge

fig, ax = plt.subplots(figsize=[10,6]);
plt.xlim(-2.5,2.5)
plt.ylim(-2.5,0.5)
ax.set_aspect(1)

# Основание, по которому скользит первое тело
base    = mlines.Line2D([-3,3],[0,0], lw = 1, color = 'k', ls = '--') 
# Подвеска тела 2
thread  = mlines.Line2D([q0[1],q0[1]+2*np.sin(q0[0])],[0,-2*np.cos(q0[0])], lw = 1, color = 'r') 
# Аннотация (время)
text    = plt.annotate('t = 0 c',[1.7,0.2])
# Тело 1 (брусок)
body1   = mlines.Line2D(q0[1]+np.array([0.2,-0.2,-0.2,0.2,0.2]),np.array([0.1,0.1,-0.1,-0.1,0.1]), lw = 2, color = 'b') 
# Тело 2 (шарик)
body2   = Wedge([q0[1]+2*np.sin(q0[0]),-2*np.cos(q0[0])], 0.1, 0, 360, fc = 'b')

def init():
    ax.add_line(base)        
    ax.add_line(thread)        
    ax.add_line(body1)    
    ax.add_patch(body2)    
    ax.add_artist(text)
    return (ax,)

def get_frame(i):        
    phi = solution.y[0, i] 
    x   = solution.y[1, i] 
    t   = solution.t[i]
    
    x1 =  x + 2*np.sin(phi)
    y1 =  0 - 2*np.cos(phi)
    
    thread.set_xdata([x, x1])
    thread.set_ydata([0, y1])
    
    body1.set_xdata(x+np.array([0.2,-0.2,-0.2,0.2,0.2]))        
    body1.set_ydata(np.array([0.1,0.1,-0.1,-0.1,0.1]))        
    
    body2.set_center([x1, y1])        
    
    text.set_text('t = {:3.2f} c'.format(t))    
    
    return (base,thread,body1,body2,text)
    
anim = FuncAnimation(fig, get_frame, init_func=init, frames = len(solution.t), interval=50, blit=False)

anim.save('animation.gif')

Вывод и интегрирование уравнений движения эллиптического мятника в Python (использование sympy, scipy)

Вывод уравнений движения

Параметры системы и обобщенные координаты

Кинематика (скорости тел)

Левая часть уравнений Лагранжа

Обобщенные силы

Приведение уравнений к форме Коши

Численное интегрирование

Графики

Анимация

ClassMech

Error

Вывод уравнений движения

Параметры системы и обобщенные координаты

Кинематика (скорости тел)

Левая часть уравнений Лагранжа

Обобщенные силы

Приведение уравнений к форме Коши

Численное интегрирование

Графики

Анимация

Templates:

Error