задача о телах

import matplotlib.pyplot as plt
from scipy import integrate
import numpy as np
 
def gravity_2(y, t, G, m1, m2):
    y1, y2, y3, y4, y5, y6, y7, y8 = y
    r3 = ((y1 - y3)**2 + (y2 - y4)**2)**(3/2)
    return [y5, y6, y7, y8, G*m2*(y3-y1)/r3, G*m2*(y4-y2)/r3, -G*m1*(y3-y1)/r3, -G*m1*(y4-y2)/r3]
 
t = np.arange(0, 5., 0.01)
y0 = [0., 0., 10., 0., 0., 3., 0., -2]
G = 1
m1 = 200
m2 = 300
 
sol = integrate.odeint(gravity_2, y0, t, args=(G, m1, m2))
plt.plot(y0[0], y0[1], '*r', y0[2], y0[3], '*g')
plt.plot(sol[:, 0], sol[:, 1])
plt.plot(sol[:, 2], sol[:, 3])
plt.title('G=1, m1=200., m2=300.')
plt.grid()
plt.show()

import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np
 
def gravity_3(y, t, G, m1, m2, m3):
    # Распаковываем входной вектор
    y1, y2, y5, y6, y3, y4, y7, y8, y9, y10, y11, y12 = y
    # Вычисляем куб расстояний
    r13 = ((y1 - y3)**2 + (y2 - y4)**2)**(3/2)
    r12 = ((y1 - y9)**2 + (y2 - y10)**2)**(3/2)
    r23 = ((y3 - y9)**2 + (y4 - y10)**2)**(3/2)
    # Возвращаем производные для системы
    return [y5, y6, 
            G*m2*(y3-y1)/r13 + G*m3*(y9-y1)/r12, 
            G*m2*(y4-y2)/r13 + G*m3*(y10-y2)/r12, 
            y7, y8, 
            -G*m1*(y3-y1)/r13 + G*m3*(y9-y3)/r23, 
            -G*m1*(y4-y2)/r13 + G*m3*(y10-y4)/r23,
            y11, y12,
            -G*m1*(y9-y1)/r12 - G*m2*(y9-y3)/r23,
            -G*m1*(y10-y2)/r12 - G*m2*(y10-y4)/r23]
 
# Начальные условия для кругового движения
# Массы и гравитационная постоянная
G = 1
m1 = m2 = m3 = 1  # Массы тел
 
R = 1  # Расстояние от центра масс до каждого тела
v = np.sqrt(G * R)  # Скорость для кругового движения
 
# Начальные условия для тел
y0 = [R, 0., 0., v, -0.5*R, np.sqrt(3)/2*R, -v, 0., -0.5*R, -np.sqrt(3)/2*R, v, 0.]
 
 
# Массив времени
t = np.linspace(0, 10, 1000)  # Временной интервал
 
# Интегрируем ОДУ
sol = odeint(gravity_3, y0, t, args=(G, m1, m2, m3))
 
# Построение графиков
plt.plot(sol[:, 0], sol[:, 1], 'r', label='Траектория m1')
plt.plot(sol[:, 4], sol[:, 5], 'g', label='Траектория m2')
plt.plot(sol[:, 8], sol[:, 9], 'b', label='Траектория m3')
plt.title('Задача трех тел с круговым движением')
plt.legend()
plt.grid()
plt.axis('equal')
plt.show()