Астрозадача №4: 

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation
from timeit import default_timer as timer

start1 = timer()
# Решение задачи многих тел
# Определяем переменную величину
seconds_in_year = 365 * 24 * 60 * 60
seconds_in_day = 60*60*24*7
years = 61
t = np.arange(0, years*seconds_in_year, seconds_in_day)

# Определяем функцию для системы диф. уравнений
def move_func(s, t):
    (x1, v_x1, y1, v_y1,
     x2, v_x2, y2, v_y2) = s

    # Динамика первого тела под влиянием второго и третьего
    dxdt1 = v_x1
    dv_xdt1 = (- G * m2 * (x1 - x2) / ((x1 - x2)**2 + (y1 - y2)**2)**1.5)
    dydt1 = v_y1
    dv_ydt1 = (- G * m2 * (y1 - y2) / ((x1 - x2)**2 + (y1 - y2)**2)**1.5)

    # Динамика второго тела под влиянием первого и третьего
    dxdt2 = v_x2
    dv_xdt2 = (- G * m1 * (x2 - x1) / ((x2 - x1)**2 + (y2 - y1)**2)**1.5)
    dydt2 = v_y2
    dv_ydt2 = (- G * m1 * (y2 - y1) / ((x2 - x1)**2 + (y2 - y1)**2)**1.5)

    return (dxdt1, dv_xdt1, dydt1, dv_ydt1,
            dxdt2, dv_xdt2, dydt2, dv_ydt2)

# Определяем начальные значения и параметры, входящие в систему диф. уравнений
G = 6.67 * 10**(-11)
Msun = 1.989*10**30
ae = 149.6*10**9
e = 0.343

x10 = 1.628*ae*(1+e)
v_x10 = 0
y10 = 0
v_y10 = 1307.03

x20 = -1.929*ae*(1+e)
v_x20 = 0
y20 = 0
v_y20 = -1548.89


s0 = (x10, v_x10, y10, v_y10,
      x20, v_x20, y20, v_y20)

m1 = 0.032*Msun
m2 = 0.027*Msun

# Решаем систему диф. уравнений
sol = odeint(move_func, s0, t)

# Строим решение в виде графика и анимируем
fig = plt.figure()
plt.style.use(['dark_background'])
fig = plt.figure(figsize=(12,12))

Q = 1.929*ae*(1+e)
x_stars = np.random.uniform(-1.3*Q, 1.3*Q, 600)
y_stars = np.random.uniform(-1.3*Q, 1.3*Q, 600)
plt.plot(x_stars/(149*10**9),y_stars/(149*10**9), 'o', ms=0.1,color= 'white')

bodys = []

for i in range(0, len(t), 10):
    body1, = plt.plot(sol[:i, 0]/(149*10**9), sol[:i, 2]/(149*10**9), '-', color='#ff6200',alpha = 0.2)
    body1_line, = plt.plot(sol[i, 0]/(149*10**9), sol[i, 2]/(149*10**9), 'o', color='#ff6200', ms = 10)

    body2, = plt.plot(sol[:i, 4]/(149*10**9), sol[:i, 6]/(149*10**9), '-', color='#ff5700',alpha = 0.2)
    body2_line, = plt.plot(sol[i, 4]/(149*10**9), sol[i, 6]/(149*10**9), 'o', color='#ff5700', ms= 12)
    
    bodys.append([body1, body1_line, body2, body2_line])

ani = ArtistAnimation(fig, bodys, interval=1)

plt.xlabel('a.e.')
plt.ylabel('a.e.')
plt.axis('equal')
# ani.save('Luhman_16.gif')
duration1 = timer() - start1
print('Время выполнения =',duration1)