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

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

The derivation of the equations of motion for an elliptical pendulum using the Python programming language and the sympy library is presented. The obtained equations are numerically integrated and an animation is created.

Рассматриваемая механическая система имеет две степени свободы: её конфигурация определяется координатой $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