Интегрирование системы дифференциальных уравнений движения системы точек в среде Python
Пример программы интегрирования уравнений движения цепи материальных точек, связанных невесомыми пружинами.
Вдоль горизонтальной прямой движутся n материальных точек с массами \(m_k\) (\(k=1,\ldots,n\)), связанные пружинами с заданными жесткостями \(c_k\) и свободными длинами \(L_k\).
Уравнения движения
Дифференциальное уравнение движение первого тела: \[m_1 \ddot{x}_1 = - c_1 (x_{1}-L_1) + c_2 (x_{2}-x_{1}-L_{2})\]
Дифференциальное уравнение движение \(k\)-го тела (\(2 \leq k \leq n-1\)): \[m_k \ddot{x}_k = - c_k (x_{k}-x_{k-1}-L_k) + c_{k+1} (x_{k+1}-x_{k}-L_{k+1}), \quad k=2,\ldots,n-1\]
Дифференциальное уравнение движение \(n\)-го тела: \[m_n \ddot{x}_n = - c_n (x_{n}-x_{n-1}-L_n)\]
Уравнения движения в форме Коши: \[\left\{ \begin{aligned} & \dot{x}_1 = V_1 \\ & \dot{V}_1 = \frac{1}{m_k} \left[- c_1 (x_{1}-L_1) + c_2 (x_{2}-x_{1}-L_{2}) \right] \\ & \dot{x}_k = V_k \quad k=2,\ldots,n-1 \\ & \dot{V}_k = \frac{1}{m_k} \left[ - c_k (x_{k}-x_{k-1}-L_k) + c_{k+1} (x_{k+1}-x_{k}-L_{k+1}), \right] \quad k=2,\ldots,n-1 \\ & \dot{x}_n = V_n \\ & \dot{V}_n = - \frac{1}{m_n} c_n (x_{n}-x_{n-1}-L_n) \\ \end{aligned} \right.\]
Программа
Подключаем библиотеки
# Массивы и матрицы
import numpy as np
# Функция численного интегрирования
from scipy.integrate import odeint
# Структура с параметрами системы
from collections import namedtuple
# Графики
import matplotlib.pylab as plt
Функция правых частей дифференциальных уравнений
def dqdt(q,t,p):
# Делим массив q (1,2n) на две части
# xv[0] = n координат x1,x2,...,xn
# xv[1] = n скоростей v1,v2,...,vn
xv = np.hsplit(q,2)
x = xv[0]
# Массив расстояний между точками
# Для первой точки - расстояние до начала координат
dx = np.diff(np.append(np.array([0]),x))
# Вычисляем силы растяжения пружин. Добавляем к этому списку справа 0
# для сохранения структуры уравнений движения последней точки,
# на которую действует сила только от одой пружины
F = np.append((dx - p.L0)*p.c,0)
# Массив сумм сил, действующих на каждую материальную точку
Fi = np.array([ -F[i]+F[i+1] for i in range(F.size-1) ])
# результат [массив скоростей, массив ускорений]
a = np.hstack((xv[1],Fi/p.m))
return a
Параметры системы
p = namedtuple("p", "n c L0 m")
# Количество тел
p.n = 5
# Массив жесткостей пружин
p.c = np.ones(p.n)*20
# Свободная длина пружин
p.L0 = 1.0
# Массы точек
p.m = np.ones(p.n)*1.0
Интегрирование уравнений движения
# Интервал интегрирования от 0 до 10 с шагом 0,01 с
ti = np.arange(0, 10.0, 0.01)
# Начальная скорость точек
v0 = np.zeros(p.n)
# Начальное положение точек
x0 = (np.random.rand(p.n)-0.5)*2*0.2+np.arange(5)+1.0
# Начальный вектор состояния
q0 = np.hstack((x0,v0))
# Запускаем процесс численного интегрирования
q = odeint(lambda q,t: dqdt(q,t,p), q0, ti)
Графики изменения координат точек
plt.figure(figsize=[11,11])
for i in range(1,6):
plt.subplot(3,2,i)
plt.plot(ti,q[:,i-1]);
plt.grid(ls=':')
plt.xlabel('t, c')
plt.ylabel('$x_{:}, м$'.format(i))
plt.tight_layout()
plt.savefig('npoints_q.png',dpi=150)
