Компьютерное моделирование движения космических тел

@hjititk · Регистрация: 27.03.2019

Студворк — интернет-сервис помощи студентам

Здравствуйте! Прошу помощи, пожалуйста. никак не могу выполнить задачу. Третье тело никак, не хочет двигаться по орбите! Буду крайне благодарен за помощь!

Задача I. Рассматривается динамика трех различных небесных тел: звезды, планеты и ее спутника. В качестве примера рассматривается Солнечная система. Масса Солнца
M1 =2⋅10^30кг. Параметры двух других тел выбираются в соответствии с индивидуальным номером варианта из ниже следующей таблицы.

Составить уравнения движения второго и третьего тела в системе отсчета, связанной с первым (самым массивным) телом. Предполагается, что движение всех тел происходит в одной плоскости.

Написать программу численного интегрирования составленных уравнений движения и построить траектории движения тел. В качестве начальных условий принять следующее: все тела находятся на одной прямой, векторы скоростей движения второго и третьего тела сонаправлены. Начальные расстояния между телами и их начальные скорости заданы в таблице выше. Там же указаны значения начальных скоростей второго и третьего тела. Исследовать отклонение орбиты планеты и спутника от круговой с течением времени, а также характер изменения их модулей скорости.

Компьютерное моделирование движения космических тел

Текущий код, при котором третье тело уходит с траектории. Так же возможно что проблема может быть в том, что должно быть гелиоцентрической системе координат. Но я не совсем понимаю это

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
 
# Гравитационная постоянная, м^3/(кг * с^2)
G = 6.67430e-11
 
# Параметры центрального тела (например, Солнце)
M1 = 2e30  # масса центрального тела, кг
 
# Параметры второго тела (например, планета)
M2 = 5.7e25        # масса второго тела, кг
R1_2 = 1.429e12    # начальное расстояние второго тела от центрального, м
V2_initial = 9.7e3 # начальная скорость второго тела, м/с
 
# Параметры третьего тела (спутник второго тела)
M3 = 1.4e23        # масса третьего тела, кг
R2_3 = 1.222e9     # начальное расстояние третьего тела от второго, м
V3_initial = 5.57e3  # начальная скорость третьего тела из таблицы, м/с
 
# Определение функции для уравнений движения
def equations(t, y):
    """
    Функция для расчета производных (ускорений и скоростей) второго и третьего тел
    под действием гравитации со стороны центрального тела и второго тела для спутника.
    """
    # Распаковка входных данных для второго тела
    x2, y2, vx2, vy2 = y[0], y[1], y[2], y[3]
    # Распаковка входных данных для третьего тела
    x3, y3, vx3, vy3 = y[4], y[5], y[6], y[7]
 
    # Расчет расстояний
    r2 = np.sqrt(x2**2 + y2**2)  # расстояние второго тела от центрального
    r3 = np.sqrt(x3**2 + y3**2)  # расстояние третьего тела (спутника) от центрального
    dx = x3 - x2                 # разница координат спутника и второго тела по x
    dy = y3 - y2                 # разница координат спутника и второго тела по y
    r3_2 = np.sqrt(dx**2 + dy**2)  # расстояние спутника от второго тела
 
    # Гравитационные ускорения второго тела относительно центрального тела
    ax2 = -G * M1 * x2 / r2**3
    ay2 = -G * M1 * y2 / r2**3
 
    # Гравитационные ускорения третьего тела относительно центрального тела и второго тела
    ax3 = -G * M1 * x3 / r3**3 - G * M2 * dx / r3_2**3
    ay3 = -G * M1 * y3 / r3**3 - G * M2 * dy / r3_2**3
 
    return [vx2, vy2, ax2, ay2, vx3, vy3, ax3, ay3]
 
# Начальные условия для второго и третьего тел
initial_conditions = [
    R1_2, 0,                      # x2, y2 - начальные координаты второго тела относительно центрального
    0, V2_initial,                # vx2, vy2 - начальные компоненты скорости второго тела
    R1_2 + R2_3, 0,               # x3, y3 - начальные координаты третьего тела относительно центрального
    0, V2_initial + V3_initial    # vx3, vy3 - начальная скорость третьего тела (спутника) вокруг второго
]
 
# Временной интервал моделирования
time_span = (0, 3.154e7 * 10)  # 10 год в секундах
time_evaluation_points = np.linspace(0, time_span[1], 5000)  # 5000 точек для траектории
 
# Численное решение задачи методом Рунге-Кутты 4-го порядка (solve_ivp)
solution = solve_ivp(
    equations,             # Функция уравнений движения
    time_span,             # Временной интервал моделирования
    initial_conditions,    # Начальные условия
    t_eval=time_evaluation_points, # Точки времени для решения
    method='RK45'          # Метод Рунге-Кутты 4-го порядка
)
 
# Построение траекторий второго и третьего тел относительно центрального тела
plt.style.use('ggplot')
 
# Построение графика с траекториями
plt.plot(solution.y[0], solution.y[1], label="Траектория второго тела", color="royalblue", linewidth=1.5)
plt.plot(solution.y[4], solution.y[5], label="Траектория спутника", color="tomato", linewidth=1.5)
plt.scatter(0, 0, color='gold', label="Центральное тело", s=100, marker='o', edgecolor='black')  # Центральное тело с контуром
 
# Настройка осей и меток с увеличенным шрифтом
plt.xlabel("x (м)", fontsize=12, fontweight='bold')
plt.ylabel("y (м)", fontsize=12, fontweight='bold')
plt.title("Траектории второго тела и спутника вокруг центрального тела", fontsize=14, fontweight='bold')
 
# Легенда за пределами графика справа
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize=10, frameon=True, fancybox=True, shadow=True, borderpad=1)
plt.grid(color='gray', linestyle='--', linewidth=0.5, alpha=0.5)
plt.axis('equal')
plt.gcf().set_size_inches(10, 6)
plt.subplots_adjust(right=0.75)
plt.show()
 
# Вывод численных данных
x2_final, y2_final = solution.y[0][-1], solution.y[1][-1]
vx2_final, vy2_final = solution.y[2][-1], solution.y[3][-1]
x3_final, y3_final = solution.y[4][-1], solution.y[5][-1]
vx3_final, vy3_final = solution.y[6][-1], solution.y[7][-1]
 
# Расчет расстояний спутника от второго тела за всё время
distances_3_2 = np.sqrt((solution.y[4] - solution.y[0])**2 + (solution.y[5] - solution.y[1])**2)
min_distance_3_2, max_distance_3_2 = np.min(distances_3_2), np.max(distances_3_2)
 
# Вывод данных
print("Конечные координаты и скорости второго тела и спутника:")
print(f"Второе тело: x = {x2_final:.2e} м, y = {y2_final:.2e} м, vx = {vx2_final:.2e} м/с, vy = {vy2_final:.2e} м/с")
print(f"Спутник: x = {x3_final:.2e} м, y = {y3_final:.2e} м, vx = {vx3_final:.2e} м/с, vy = {vy3_final:.2e} м/с")
print()
print("Минимальные и максимальные расстояния спутника от второго тела:")
print(f"Спутник: min = {min_distance_3_2:.2e} м, max = {max_distance_3_2:.2e} м")

@hjititk · 06.11.2024, 19:08 **[ТС]**

Я так понимаю вся загвоздка кроется здесь?

Python
1
2
    ax3 = -G * M1 * x3 / r3**3 - G * M2 * dx / r3_2**3
    ay3 = -G * M1 * y3 / r3**3 - G * M2 * dy / r3_2**3

Python
1
    R1_2 + R2_3, 0, V2_initial + V3_initial * np.cos(np.pi / 2), V3_initial * np.sin(np.pi / 2)

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
 
# Гравитационная постоянная, м^3/(кг * с^2)
G = 6.67430e-11
 
# Параметры центрального тела Солнце
M1 = 2e30  # масса центрального тела, кг
 
# Параметры второго тела (планета)
M2 = 5.7e25        # масса второго тела, кг
R1_2 = 1.429e12    # начальное расстояние второго тела от центрального, м
V2_initial = 9.7e3 # начальная орбитальная скорость второго тела вокруг центрального, м/с
 
# Параметры третьего тела (спутник второго тела)
M3 = 1.4e23        # масса третьего тела, кг
R2_3 = 1.222e9     # начальное расстояние третьего тела от второго, м
V3_initial = 5.57e3  # начальная скорость третьего тела относительно второго, м/с
 
# Определение функции для уравнений движения
def equations(t, y):
    x2, y2, vx2, vy2 = y[0], y[1], y[2], y[3]
    x3, y3, vx3, vy3 = y[4], y[5], y[6], y[7]
 
    r2 = np.sqrt(x2**2 + y2**2)
    r3 = np.sqrt(x3**2 + y3**2)
    dx = x3 - x2
    dy = y3 - y2
    r3_2 = np.sqrt(dx**2 + dy**2)
 
    ax2 = -G * M1 * x2 / r2**3
    ay2 = -G * M1 * y2 / r2**3
 
    ax3 = -G * M1 * x3 / r3**3 - G * M2 * dx / r3_2**3
    ay3 = -G * M1 * y3 / r3**3 - G * M2 * dy / r3_2**3
 
    return [vx2, vy2, ax2, ay2, vx3, vy3, ax3, ay3]
 
# Начальные условия для второго и третьего тел
initial_conditions = [
    R1_2, 0, 0, V2_initial,
    R1_2 + R2_3, 0, V2_initial + V3_initial * np.cos(np.pi / 2), V3_initial * np.sin(np.pi / 2)
]
 
# Временной интервал моделирования
time_span = (0, 3.154e7 * 10)
time_evaluation_points = np.linspace(0, time_span[1], 10000)
 
# Численное решение задачи методом Рунге-Кутты 4-го порядка (solve_ivp)
solution = solve_ivp(equations, time_span, initial_conditions, t_eval=time_evaluation_points, method='RK45')
 
# Построение траекторий второго и третьего тел относительно центрального тела
plt.style.use('ggplot')
plt.plot(solution.y[0], solution.y[1], label="Траектория второго тела", color="royalblue", linewidth=1.5)
plt.plot(solution.y[4], solution.y[5], label="Траектория спутника", color="tomato", linewidth=1.5)
plt.scatter(0, 0, color='gold', label="Центральное тело", s=100, marker='o', edgecolor='black')
 
plt.xlabel("x (м)", fontsize=12, fontweight='bold')
plt.ylabel("y (м)", fontsize=12, fontweight='bold')
plt.title("Траектория второго тела и спутника вокруг центрального тела", fontsize=14, fontweight='bold')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize=10, frameon=True, fancybox=True, shadow=True, borderpad=1)
plt.grid(color='gray', linestyle='--', linewidth=0.5, alpha=0.5)
plt.axis('equal')
plt.gcf().set_size_inches(10, 6)
plt.subplots_adjust(right=0.75)
plt.show()

@Aael · 08.11.2024, 18:32

hjititk,

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from math import sqrt
 
def gravity(x1, y1, m1, x2, y2, m2):
    dx, dy = x2 - x1, y2 - y1
    r = sqrt(dx * dx + dy * dy)
    f = G * m2 / r / r / r
    # f = G * (m1 + m2) / r / r / r
    return f * dx, f * dy
 
def clamp_abs(*seq, limit=1e6):
    '''на случай если планета улетит в бесконечность'''
    return [min(max(-limit, x), limit) for x in seq]
 
Gm = 1 # гигаметр (млн. км)
km = Gm / 1000000
m = km / 1000
 
kg = 1e-23
 
year = 1
day = year / 365
hour = day / 24
sec = hour / 60 / 60
 
G = 4405039 / 660000 * 1e-11 * m**3 / kg / sec**2
 
# Солнце
M1 = 2e30 * kg
# Сатурн
M2 = 5.7e26 * kg
R12 = 1429e6 * km
V2 = 9.7 * km / sec
T21 = 29.5 * year # период обращения Сатурна вокруг Солнца
# Титан
M3 = 1.4e23 * kg
R23 = 1222e3 * km
V3 = 5.57 * km / sec
T23 = 16 * day # период обращения Титана вокруг Сатурна
 
def fn(t, y):
    x2, y2, vx2, vy2 = y[:4]
    x3, y3, vx3, vy3 = y[4:]
    ax2, ay2 = gravity(x2, y2, M3, 0, 0, M1)
    ax31, ay31 = gravity(x3, y3, M3, 0, 0, M1)
    ax32, ay32 = gravity(x3, y3, M3, x2, y2, M2)
    return clamp_abs(vx2, vy2, ax2, ay2, vx3, vy3, ax31 + ax32, ay31 + ay32)
 
fn0 = [
    R12, 0, 0, V2,
    R12 + R23, 0, 0, V2 + V3,
]
 
resolution = 5000
plt.style.use('ggplot')
fig, ((p1, p2, p3)) = plt.subplots(1, 3)
fig.suptitle('Траектории движения Сатурна и Титана')
 
dt = (0, 1 * T21)
dts = np.linspace(0, dt[1], resolution)
s_large = solve_ivp(fn, dt, fn0, 'RK45', dts, max_step = day/2)
print(s_large.message)
 
p1.set_title('За период Сатурна\n(траектории сливаются,\nможно приблизить)', fontsize=10)
line1, = p1.plot(s_large.y[0], s_large.y[1], label="Сатурн", color="royalblue", linewidth=1.5)
line2, = p1.plot(s_large.y[4], s_large.y[5], label="Титан", color="tomato", linewidth=1.5)
line3  = p1.scatter(0, 0, color='gold', label="Солнце", s=100, marker='o', edgecolor='black')
p1.set_xlabel("x (млн.км)")
p1.set_ylabel("y (млн.км)")
p1.axis('equal')
 
times = 10
dt = (0, times * T23)
dts = np.linspace(0, dt[1], resolution)
s_small = solve_ivp(fn, dt, fn0, 'RK45', dts, max_step = hour)
print(s_small.message)
 
p2.set_title(f'За {times} периодов Титана\nцентр ск - Солнце', fontsize=10)
p2.plot(s_small.y[0], s_small.y[1], label="Сатурн", color="royalblue", linewidth=1.5)
p2.plot(s_small.y[4], s_small.y[5], label="Титан", color="tomato", linewidth=1.5)
p2.set_xlabel("x (млн.км)")
p2.axis('equal')
 
p3.set_title(f'За {times} периодов Титана\nцентр ск - Сатурн', fontsize=10)
p3.scatter(0, 0, color='royalblue', label="Сатурн", s=100, marker='o', edgecolor='black')
p3.plot(s_small.y[4] - s_small.y[0], s_small.y[5] - s_small.y[1], label="Титан", color="tomato", linewidth=1.5)
p3.set_xlabel("x (млн.км)")
p3.axis('equal')
 
fig.legend(handles=[line1, line2, line3], loc='lower center',ncol=3)
plt.subplots_adjust(top=0.8, bottom=0.2)
plt.show()

@hjititk · 08.11.2024, 18:42 **[ТС]**

Большое Вам спасибо!

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Web (WebAssembly): Обработчик клика мыши в браузере ПК и касания экрана в браузере на мобильном устройстве 8Observer8 02.02.2026 Содержание блога Для начала пошагово создадим рабочий пример для подготовки к экспериментам в браузере ПК и в браузере мобильного устройства. Потом напишем обработчик клика мыши и обработчик. . .	Философия технологии iceja 01.02.2026 На мой взгляд у человека в технических проектах остается роль генерального директора. Все остальное нейронки делают уже лучше человека. Они не могут нести предпринимательские риски, не могут. . .	SDL3 для Web (WebAssembly): Вывод текста со шрифтом TTF с помощью SDL3_ttf 8Observer8 01.02.2026 Содержание блога В этой пошаговой инструкции создадим с нуля веб-приложение, которое выводит текст в окне браузера. Запустим на Android на локальном сервере. Загрузим Release на бесплатный. . .	SDL3 для Web (WebAssembly): Сборка C/C++ проекта из консоли 8Observer8 30.01.2026 Содержание блога Если вы откроете примеры для начинающих на официальном репозитории SDL3 в папке: examples, то вы увидите, что все примеры используют следующие четыре обязательные функции, а. . .
SDL3 для Web (WebAssembly): Установка Emscripten SDK (emsdk) и CMake для сборки C и C++ приложений в Wasm 8Observer8 30.01.2026 Содержание блога Для того чтобы скачать Emscripten SDK (emsdk) необходимо сначало скачать и уставить Git: Install for Windows. Следуйте стандартной процедуре установки Git через установщик. . . .	SDL3 для Android: Подключение Box2D v3, физика и отрисовка коллайдеров 8Observer8 29.01.2026 Содержание блога Box2D - это библиотека для 2D физики для анимаций и игр. С её помощью можно определять были ли коллизии между конкретными объектами. Версия v3 была полностью переписана на Си, в. . .	Инструменты COM: Сохранение данный из VARIANT в файл и загрузка из файла в VARIANT bedvit 28.01.2026 Сохранение базовых типов COM и массивов (одномерных или двухмерных) любой вложенности (деревья) в файл, с возможностью выбора алгоритмов сжатия и шифрования. Часть библиотеки BedvitCOM Использованы. . .	SDL3 для Android: Загрузка PNG с альфа-каналом с помощью SDL_LoadPNG (без SDL3_image) 8Observer8 28.01.2026 Содержание блога SDL3 имеет собственные средства для загрузки и отображения PNG-файлов с альфа-каналом и базовой работы с ними. В этой инструкции используется функция SDL_LoadPNG(), которая. . .

@hjititk 0 / 0 / 0 Регистрация: 27.03.2019 Сообщений: 17
	08.11.2024, 18:42 [ТС]
	Большое Вам спасибо! 0

Компьютерное моделирование движения космических тел

Решение