Движение пары взаимодействующих частиц

Две частицы притягиваются с силой пропорциональной \(1/r^2\). Построим траектории их движения. Физическая запись задачи содержит два векторных ДУ уравнения второго порядка:

\[ \begin{align}\begin{aligned}m_1 \frac{d^2\vec{r_1}}{dt^2}=\vec{F}_1,\\m_2 \frac{d^2\vec{r_2}}{dt^2}=\vec{F}_2.\end{aligned}\end{align} \]

Для формализации задачи и записи ее в форме \(\frac{d}{dt}\vec{x}(t)=\vec{f}(\vec{x},t)\) необходимо определить вектор формальных переменных: \(\vec{x}=[t, x_1, y_1, v_{1x},v_{1y}, x_2, y_2, v_{2x}, v_{2y}]\). Для удобства время считается одной из переменных.

Производные соответствующих переменных дают вектор правых частей уравнения: \(\vec{f}=\frac{d\vec{x}}{dt}=[1, v_{1x}, v_{1y}, F_{1x}/m_1, F_{1y}/m_1, v_{2x}, v_{2y}, F_{2x}/m_2, F_{2y}/m_2]\). Задаем функцию, вычисляющую вектор производных (правая часть ДУ):

def force(x):
#  vector of data       t, x1, y1,vx1,vy1, x2, y2,vx2,vy2
#  vector of derivative 1,vx1,vy1,fx1,fy1,vx2,vy2,fx2,fy2

 r=((x[1]-x[5])**2+(x[2]-x[6])**2)**0.5
 fx1=-(x[1]-x[5])/r**2
 fy1=-(x[2]-x[6])/r**2
 fx2=-fx1
 fy2=-fy1
 f=np.array([1., x[3], x[4], fx1,fy1, x[7], x[8], fx2,fy2])
 return f

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

def step(f, x, h): #
  kt1=f(x)*h
  kt2=f(x+kt1/2.)*h
  kt3=f(x+kt2/2.)*h
  kt4=f(x+kt3)*h
  xn=x+(kt1+2.*kt2+2.*kt3+kt4)/6. 
  return xn

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

####### main part #######################
#  vector of data t,x1,y1,vx1,vy1,x2,y2,vx2,vy2
x=np.array([0, 1., 0, 0, 1., -1., 0, 0, -1.]) 
h=0.1

py.ion()
py.xlabel('X')
py.ylabel('Y')

for nt in range(1000):
  x=step(force, x, h)
  py.plot(x[1], x[2], 'go') #
  py.plot(x[5], x[6], 'bo') #
  py.draw()

Выполнить описанный пример