Пользователь
0,0
рейтинг
16 октября 2013 в 12:34

Разработка → Создаем симулятор солнечной системы из песочницы

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

Волею судьбы на глаза попался замечательный цикл статей о создании игры-платформера на Питоне
тут и тут.
Я решил взяться за один старый проект. За симулятор движения тел под действием сил гравитации.

Что из этого вышло читайте дальше.

Часть первая. Теоритическая

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

Введем систему координат.

Пускай наша система состоит из двух тел:
1. массивной звезды массой М и центром (x0, y0)
2. легкой планеты массой m, с центром в точке (x, y), скоростью v = (vx, vy) и ускорением a = (ax, ay).

Когда нам удастся разобрать этот случай, мы легко перейдем к сложным системам со взаимным влиянием звезд и планет друг на друга. Сейчас же речь пойдет о самом простом.

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

ax = G * M * (x0-x) / r^3
ay = G * M * (y0-y) / r^3

Это позволяет составить алгоритм перемещения планеты в поле гравитации звезды:

1. Перед началом задаем начальное положение планеты (x, y) и начальную скорость (vx, vy)
2. На каждом шаге вычисляем новое ускорение по формуле выше, после этого пересчитываем скорость и координаты:

vx := vx + T * ax
vy := vy + T * ax

x := x + T * vx
y := y + T * yx

Осталось разобраться с константами G и T. Положим G = 1. Для нашей задачи это не так важно. Параметр T влияет на точность и скорость вычислений. Тоже положим 1 для начала.

Часть вторая. Практическая

Итак, моя первая программа на Питоне. При этом еще раз хочется поблагодарить Velese за практическое руководство.

import pygame, math
from pygame import *
from math import *

WIN_WIDTH = 800
WIN_HEIGHT = 640
PLANET_WIDTH = 20
PLANET_HEIGHT = 20
DISPLAY = (WIN_WIDTH, WIN_HEIGHT)
SPACE_COLOR = "#000022"
SUN_COLOR = "yellow"
PLANET_COLOR = "blue"

#Sun position
X0 = WIN_WIDTH // 2
Y0 = WIN_HEIGHT // 2
#Sun mass
M0 = 5000
#Stop conditions
CRASH_DIST = 10
OUT_DIST = 1000

def main():
    #PyGame init
    pygame.init()
    screen = pygame.display.set_mode(DISPLAY)
    pygame.display.set_caption("Solar Mechanics v0.1")
    
    #Space init
    bg = Surface((WIN_WIDTH,WIN_HEIGHT))
    bg.fill(Color(SPACE_COLOR))     
    draw.circle (bg, Color(SUN_COLOR), (X0, Y0), 10)
                    
    #Timer init                     
    timer = pygame.time.Clock()
    
    #Planet init
    planet = Surface((PLANET_WIDTH, PLANET_HEIGHT))
    planet.fill(Color(SPACE_COLOR))
    draw.circle (planet,
        Color(PLANET_COLOR),
        (PLANET_WIDTH // 2, PLANET_HEIGHT // 2),
        5)

    #Planet to Sun distance
    r = 0.0
    #Initial planet pos, speed and accel
    x = 100.0
    y = 290.0
    vx = 0.1
    vy = 1.5
    ax = 0.0
    ay = 0.0

    done = False
    while not done:
        timer.tick(50)
        for e in pygame.event.get():
            if e.type == QUIT:
                done = True
                break        

        r = sqrt((x - X0)**2 + (y - Y0)**2)
        
        ax = M0 * (X0 - x) / r**3
        ay = M0 * (Y0 - y) / r**3

        #New spped based on accel
        vx += ax
        vy += ay

        #New pos based on speed
        x += vx
        y += vy

        screen.blit(bg, (0, 0))      
        screen.blit(planet, (int(x), int(y)))
        pygame.display.update()     

        if r < CRASH_DIST:
            done = True
            print("Crashed")
            break
        if r > OUT_DIST:
            done = True
            print("Out of system")
            break

    #Farewell
    print (":-)")

if __name__ == "__main__":
    main()


Так выглядит наша система после некоторого времени симуляции



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

Сейчас я занимаюсь поиском интересных сценариев системы и небольшими улучшениями интерфейса.

Вот пример того, что на данный момент в разработке:


Если эта заметка встретит положительные отзывы, обещаю продолжить рассказ о более новой версии.

Update:

1. Я благодарен всем комментаторам за критические замечания. Они дают большую пищу для размышлений.

2. Проект вырос. Все тела уже независимы, влияют друг на труга в соответствии с законом всемирного тяготения.Подсчитывается N^2 взаиможействий.
Сейчас есть возможность хранить конфигурации звездной системы во внешних файлах и выбирать на старте
Код тут
Запускать так: python3.3 main.py -f <имя конфигурации>.ini
Различные конфигурации — там же.

3. Благодаря комментариям удалось найти и устранить главную недоработку — метод вычисления координат.
Сейчас используется метод Рунге-Кутты. По мере прочтения «Нежестких задач» буду осваивть новые методы.
@RedTroll
карма
9,0
рейтинг 0,0
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Спецпроект

Самое читаемое Разработка

Комментарии (69)

  • +2
    Вот так вот и рождаются игры типа Kerbal Space Program, Free Orion, Solar (2).
    • +1
      в Kerbal Space Program не учитывается одновременное влияние двух тел, следовательно, невозможны точки Лагранжа. К слову сказать, задача о движении n тел под воздействием гравитации невероятно сложна для аналитического решения уже при n=3.
      • 0
        невероятно сложна для аналитического решения уже при n=3

        Вы хотели сказать, в общем случае аналитически неразрешима?
        • 0
          Не позорьтесь, о Зундмане даже на Вики написано.
    • 0
      Тогда уж, скорее, Universe Sandbox.
  • +8
    Чего нельзя изобразить
    books.google.ru/books?id=i7V2LXYdC-0C&pg=PA145
    • +1
      Залип минут на 20
  • 0
    То же такие мысли были. Я правда сразу за 3D взялся, время от времени кодю. Но у меня ещё ничего не летает. А вам удачи в проекте)
    • 0
      Спасибо. Возможно, я дойду до 3D. Пока это для меня сложно.

      Удачи и Вам
    • +1
      Фишка в том, что задача даже в 2D настолько сложна, а процесс настолько красив, что можно смело забыть про третье измерение и считать, что божественным произволом все тела оказались в одной плоскости, z=0.
      • 0
        Если мы говорим о симуляторе солнечной системы, то без третьего измерения в конце концов не обойтись.
        • 0
          Если учитывать только планеты — то третье измерение лишнее…
  • 0
    Вам мой положительный отзыв: развивайте проект! Вы мне напомнили давнюю мою мечту, о которой я уже давно забыл, а сейчас она возродилась, так что буду следить.

    А можно уже новую версию опробовать?
    • 0
      Спасибо за отзыв. Новая версия тут:

      github.com/variostudio/spacesim

      Исполняемый файл: main.py
      Есть ключ -h
      Есть ключ -f для загрузки .ini файла с конфигурацией системы
      Есть несколько конфигураций. Пробуйте. Пишите что получилось.
      Возможно, что-то криво закоммичено. Пишите, вечером поправлю
      • 0
        Не получилось поставить ConfigParser почему–то (Win7/Python2.6), но развил версию приведёную здесь до двух планет и звезда тоже обращается вокруг общего центра масс. Захватывает… Глядишь, разовьётся в вот такое:
        • 0
          Досадно…
          Я, как начинающий в Питоне, писал сразу под версию 3.3
          Может, в этом проблема?

          Теоритически вы можете заменить ConfigParser на что-то из набора Python2.6

          Там есть несколько готовых систем.
          • 0
            Не очень понятно – тоже новичёк.
            Кстати, чтобы поделится программой с теми кто не имеет Питона и сидит на виндах, можете скомпиливать в exe с помощью py2exe, а лучше с pywin32 (но может это только из–под винды будет работать...)

            П.С. Планет уже 2, полная нестабильность, пока одну не вышвырнет из системы… или все разлетаются. Такие пируэты выдают! Как вы и писали, вся фишка в начальных данных, можно поймать очень красивую комбинацию.
            • 0
              Я попробую откомпилировать вечером в.ехе
              Да, иногда бывает сложно поймать стабильную комбинацию :-) Чем больше планет тоем тоньше грань.

              Попробуйте Питон 3.3, если есть возможность.
              Все должно запуститься.
              • 0
                Ага. Прочитал, что до 3.3 configparser назывался ConfigParser

                The ConfigParser module has been renamed to configparser in Python 3. The 2to3 tool will automatically adapt imports when converting your sources to Python 3

                Может, это поможет.
                • 0
                  Получилось с Питоном версии 3.2, версия 3.3 не вязалась с pygame. Спасибо!
                  • 0
                    Удалось разные конфигурации запустить (ключ -f и имя .ini файла)?
                    • 0
                      Да, всё отлично работает!

                      Кстати, в своей версии я выводил ещё текстовую информацию о расстояниях в окно. Возможно, кому-то будет интересно наблюдать как меняются параметры орбиты той планеты, на которую щёлкнули мышью. Также, интересно поэкспериментировать с тяжёлыми планетами, но тогда нужен режим (с возможностью отключить) где звезда не стационарна, а тоже динамична как и планеты. Я уже написал, что в своей версии я запускал планеты с массами близкими к звёздной — красота :)
  • 0
    Видео есть? Если оно летает уже
    • 0
      Летает.
      Как сделать видео со скрина я пока не знаю. KUbuntu.
      Но в скором времени займусь.
  • 0
    Делал, что-то подобное на js + canvas. Планеты у меня можно добавлять кликом и потянуть для задания размера.
    • 0
      А как задавались начальные скорости для планет?
      • 0
        Рандомом, сейчас чуток переделал на ваши формулы, чуток изменю логику.
        • 0
          В начальных скоростях вся соль. Все великие столкновения планет или побеги из системы — суть начальных скоростей. Рандомом нельзя…

          Насчет логики. Сначала нужно пересчитать ВСЕ ускорения а ТОЛЬКО ПОТОМ — все скорости и координаты.
          Нельзя считать одновременно ускорения, скорости и координаты в цикле по объктам, иначе новые значения координат будут влиять на ускорения других объектов.
          Об этом будет вторая статья.
          • 0
            спасибо, а то и правда сумашествие какое-то происходит
          • 0
            Насчет логики. Сначала нужно пересчитать ВСЕ ускорения а ТОЛЬКО ПОТОМ — все скорости и координаты.
            Нельзя считать одновременно ускорения, скорости и координаты в цикле по объктам, иначе новые значения координат будут влиять на ускорения других объектов.

            Ну вы вот сами запутались.

            Законы Ньютона для N тел напрямую дают систему из 3N дифференциальных уравнений второго порядка вида:
            x_i'' = f_i(x_1, .. ,x_n, x_1', .., x_n', t)
            Здесь разные x — это какая-то из трёх координат каждого тела, всего 3N координат — итого 3N таких уравнений. x — координаты, x' — скорости, x'' — ускорения (ну, т.е. скорости и ускорения — первые и вторые производные координат)/
            Функции в правой части на самом деле не зависят от скоростей и времени, но я просто хотел записать общий случай. Они, собственно говоря, сумма гравитационных взаимодействий, зависят только от координат.

            Каждое такое уравнение эквивалентно двум уравнениям первого порядка, а вся система эквивалентна системе из 6N уравнений первого порядка вида (ещё я учёл, что гравитация от скорости и времени не зависит):
            x_i' = v_i
            v_i' = f(x_1, .., x_n)
            (т.е. 3N уравнений первого типа и 3N второго)

            здесь x — это координаты, x' = v — скорости, v' — ускорения, а в целом если задать ещё начальные условия (3N координат и 3N скоростей), получим типичную задачу Коши.

            Если производные в левой части заменить конечными разностями, мгновенно получаем метод Эйлера: x' = (x{k+1}-x{k})/dt, отсюда x{k+1} = x{k}+dt*f(...) (6N таких уравнений — по одному для каждого x и v) Значит, теперь можно работать методами Рунге-Кутта и прозводными, Адамса и производными — вам осталось построить явный вид функций правой части.

            В любом случае, для явных методов семейства Рунге-Кутта вы будете иметь выражение вида x{k+1} = R[x{k}, dt, f(...)] (выражение метода Эйлера выше — частный случай). Это выражение будет весьма сложным, зато оно универсальное для любых задач интегрирования (и уже миллион раз реализовано в разных библиотеках). Некоторые методы ещё будут модифицировать шаг и давать помимо очередного значения погрешность.

            Что нужно понять из этого комментария.
            1. Скорости и координаты с точки зрения математики здесь практически ничем друг от друга не отличаются. То, что уравнений всего 6N — прямое следствие того, что система N тел имеет фазовое пространство размерности 6N, а эти размерности могут быть, например, декартовыми координатами и скоростями.
            2. Системы ДУ решаются как одно ДУ, под интегралом и в аргументе правой части задачи Коши которого стоит вектор. Этот вектор — и есть штука, которую вы пытаетесь посчитать.
            3. Ну, как положено, считая следующее значение этого вектора, нужно хранить целиком предыдущее значение. Если вспомнить, что вектор имеет вид {3N координат, 3N скоростей}, то вы храните предыдущие координаты и скорости. Хранить и пересчитывать ускорения вообще не нужно — они после дискретизации задачи (перехода к конечным разностям) вообще пропадают из всех формул. (Если строго говорить, их в формулах и не было — только в виде второй производной от координаты или первой от скорости.)
            • 0
              Да. Я именно про п.3 говорил.
              Попервах у меня наслаивалось следующее значение на предыдущее. Простой баг был :-)
            • 0
              Запишу задачу так:
              X = \begin{bmatrix} x_1 \\ .. \\ x_n \\ v_1 \\ .. \\ v_n \end{bmatrix}
              X — это тот самый вектор
              X' = F(X, t)
              А это — уравнение, которое мы решаем. Для двух тел правая часть имеет вид (от t она явно не зависит):
              F(X) = \begin{bmatrix} k m_2 (x_1-x4) \\ k m_2 (x_2-x_5) \\ k m_2 (x_3-x_6) \\ k m_1 (x_4-x_1) \\ k m_1 (x_5-x_2) \\ k m_1 (x_6-x_3) \\ x_1 \\ .. \\ x_6 \end{bmatrix}
              где k — это общий знаменатель
              k=\((x_1-x_4)^2+(x_2-x_5)^2+(x_3-x_6)^2\)^{-\frac{3}{2}}
              Мог напутать разве что со знаками :) И забыл гравитационную постоянную, конечно же. Входит в k в виде ещё одного сомножителя.
              x1..x3 — координаты первого тела, x4..x6 — второго тела, x7..x9 — скорости первого, x10..x12 — второго.
      • 0
        Если вам интересно, то могу написать пост, как смоделировать закон Кулона (ну или всемирного тяготения — просто константы другие :) ). Я когда-то писал на Хабре один комментарий на эту тему, но народ потом переключился на обсуждение физики, решение диф. уравнений и так далее.
        Посмотрите, там есть две старые видюшки моего JavaScript + Canvas решения. Скорей всего то, что d4rkr00t делал давно.
        • 0
          вообще да, было бы очень интересно, у меня то там совсем все просто было codepen.io/d4rkr00t/pen/hItgE движение по окружности + что-то типа силы грацитации.
        • 0
          Спасибо за комментарий и ссылку.

          Да, очень интересно. Был бы рад увидеть Вашу статью.
          И на программу d4rkr00t тоже интересно взглянуть :-)
          • 0
            старая версию по ссылке выше, на codepen.io
  • 0
    Лет 5 назад пользовался подобной программой, в ней тела при столкновении слипались, и можно было симулировать образование планет. Неплохо бы подобную вещь, но со сбалансированным слипанием, зависящим от притяжения.
    • 0
      Образование новых планет — интересная мысль. Я уже думал про слипание, но немного пугает, что при лобовом ударе может происходить раскалывание.

      В любом случае запишу себе в планы.

      Спасибо.
  • +1
    Добрый день. Несколько замечаний по алгоритму.
    Закон всемироного тяготения у вас записан не совсем верно: во-первых, это сила, во-вторых у нее две массы (массы обоих тел). Конечно, если силу подставить во второй закон Ньютона и расписать на оси, так и получится, но формулировка у вас дана неверная.
    Отсюда вывод — сделайте его как силу, а в уравнения движения (которые в вашем случае являются интегрированием методом Эйлера, который имеет только одним преимуществом — простотой, но недостатков у него вагон) подставить силу, деленную на собственную массу. Это позволит улучшить модель и сделать взаимное влияние планет друг на друга (а оно очень большое), да и можно будет просто добавить, например, Луну или другие спутники. В идеале это еще позволит учитывать кучу возмущений, но вам вряд ли пока это нужно.
    Второе — познакомьтесь с кеплеровскими элементами орбиты. Опять две причины: а — они позволяют не париться о выборе начальных условий (собственно, вы задаете элементы орбиты, а дальше они пересчитываются в XYZ, vx, vy, vz), б — в небесной механики координаты тел задаются как раз ими, так что если захотите брать реальные цифры — обязательно с ними столкнетесь.
    • 0
      Спасибо за комментарий.

      Сейчас я пишу вторую часть статьи, где учтено взаимное влияние тел друг на друга. В примерах есть симуляция системы «Звезда — планета — спутник». Даже в двух вариантах.

      Можно ли будет обратиться к Вам для для комментария обновленной версии?

      Кстати, вопрос. В новой версии я пользуюсь этим же алгоритмом, но вычисляемом для каждой пары тел в обе стороны.
      Т.е. есть тела 1...N
      Я вычисляю ускорение при действии тела №i на тело №j, потом наоборот, тела №j на тело №i
      Все суммируется. Получается N^2 вычислений для системы из N тел. Каждое взаимодействие учтено.
      Этот подход имеет право на существование?

      Еще раз спасибо за конструктивный комментарий.
      • 0
        Не очень понял с ускорениями.
        Я бы советовал двигаться в сторону второго закона Ньютона. А именно — ускорение считать как Ai = Sum(Fij)/m, где Fij — сила взаимного притяжения тел (Fij = G * mi * mj * (Xi — Xj) / (rij ^ 3), тут записано в векторной форме).
        Соответственно Vi = Vi + dt * Ai, Xi = Xi + dt * Vi; где dt — шаг интегрирования.
        Эм. Путаница с индексами, так что считайте, что это присвоение, а не равенство, чтобы не рассписывать V_i, t каждый раз.
        При интегрировании в реальном времени dt равен времени с прошлого рассчета (если у вас 24 кадра в секунду, то оно равно 1/24, но лучше засекать время с предыдущей итерации).

        Да, ко мне можно обращаться, всегда с радостью помогу.
        • 0
          Не очень понял с ускорениями.
          Я бы советовал двигаться в сторону второго закона Ньютона. А именно — ускорение считать как Ai = Sum(Fij)/m, где Fij — сила взаимного притяжения тел (Fij = G * mi * mj * (Xi — Xj) / (rij ^ 3), тут записано в векторной форме).

          Я ж точно так и делал. Только сразу сократил m и оперировал не понятием «сила», а понятием «ускорение» :-)

          Мне уже подсказали воспользоваться другими численными методами. Возможно, я так и поступлю, если удастся осилить. Университетский курс численных методов требует много времени. Школьный курс физики намного менее требовательный.
      • 0
        Вообще хорошо бы вам почитать первый том Ландау-Лифшица. Сложно, да. Но там научат строить строго систему уравнений движения, дифференциальных, в данном случае они будут второго порядка. Ну да из системы N второго порядка сделать 2N первого — раз плюнуть. Это будет универсальный метод для любых количеств любых тел.

        А эти уравнения можно интегрировать прямо сразу численно.
  • +1
    Использованный метод численного дифференцирования фактически представляет собой метод Эйлера первого порядка, и, в случае эллиптических орбит, будет быстро набегать систематическая ошибка — планеты будут потихоньку съезжать «от Солнца».
    Традиционно в программах такого рода используют метод Рунге-Кутты 4-го порядка. Он лишь чуть сложнее метода Эйлера в программировании, но намного точнее.
    • 0
      Точно, так и есть…

      Буду осваивать новый метод. Есть над чем поработать теперь.

      Огромное спасибо.

      • +1
        Дополню. Во-первых, используется не сам метод Ронге-Кутта, а его модификация — Дормана-Принса. Он хорош тем, что пересчитывает шаг и где можно шагать без потерь огромными шагами так и делает, а где нужно точно проинтегрировать — уменьшает шаг. Во-вторых, очень удобно в подобных, да и во всех вообще случаях. Но если вам вдруг взбредет симулировать движение планет на продолжительный период времени, стоит обратить внимание на орбито-устойчивые методы. Они позволяют интегрировать сотни, тысячи, миллионы оборотов с достаточной точностью. С ходу ссылку найти не удалось, но могу подсказать литературу, если нужно.
        • 0
          ДП — няшечка. Обычно гоняют 5(4), но есть он и более высокого порядка — 8(7)
          • 0
            У меня никогда не хватало терпения вбить всю таблицу для 8(7) порядка))
            Я вот давно хочу еще попробовать метод Адамса, но руки не доходят. Сейчас буду комплексировать GPS/ГЛОНАСС и радио-баровысотомер (модели, конечно), попробую написать Адамса или еще что-то из экзотики)
            • 0
              +1 про ДП8(7) :) я по-моему его вбил, считало неправильно, искать ошибку… хрен найдёшь её.

              Есть ещё неявные методы семейства РК, там свои плюсы (устойчивее) и минусы (реализовывать посложней будет).

              Адамс хорош, когда реально должно получаться что-то очень близкое к многочленам. Он же аппроксимирует многочленами. Если точно решение — многочлен порядка K, то метод Адамса любоого порядка P>K должен давать 0 в старших членах, т.к. он точно найдёт K+1 коэффициентов этого многочлена. Отсюда пляшем выбор метода интегрирования.
    • +1
      За РК четвёртого — минус вам. Ну хватит уже его насиловать. Морально устарел РК 4-го порядка, никто его не использует, только некоторые, кто «что-то слышал».
  • 0
    Когда нам удастся разобрать этот случай, мы легко перейдем к сложным системам со взаимным влиянием звезд и планет друг на друга.


    На этой фразе наверное все физики рассмеялись.
    • 0
      В рамках методов прямого численного дифференцирования N-body систем, как у автора — действительно легко.
      Вот если аналитически — тут да, уже решение задачи трёх тел весьма проблематично.
      • 0
        Да нууууу? Легко, да?

        Как вы думаете, в системе из 10 тел, далеко вы уйдёте методом Эйлера, чтоб погрешности не выросли хотя бы на 10%?

        Очень недалеко.
        • 0
          Даже Эйлером, на пару сотен характеристических времен пересечения — вполне себе. Рецепт прост — следите за общими интегральными характеристиками системы, такими как интеграл площадей. Видите, что ползёт в сторону — режьте шаг :)
          • 0
            Тупо x''=-x уже не интегрируется Эйлером. Уже после нескольких колебаний видно, что амплитуда растёт — а должен был получиться синус. Ну можно конечно обрезать шаг, но сколько это будет считаться?

            И это — всего одна переменная.
  • +1
    > Чтобы решить задачу, нужно сначала четко себе ее представить.
    Вот-вот.

    Рекомендую ознакомится с исходниками давно построенного и проверенного велосипеда под названием "Симулятор Солнечной системы".
    А именно, что б долго не искать, с классом Orbit и, например, его использованием в методе Body::getPosition().

    p.s. Ваша заметка плохая. Теории нет, решения заявленной задачи нет. То, что вы накодили — кривое и, строго говоря, нерабочее… Вообще оставляет впечатление «нам вчера в классе рассказали про питон и зацените, что я сегодня набросал».
    Не надо такую же вторую писать.
    • 0
      Да, тут я вынужден согласиться. На момент написания я помнил только закон всемирного тяготения и 2-й закон Ньютона. Пробел в этой области
      Но теперь я знаю что есть метод Эйлера и метод Рунге-Кутты.

      Спасибо за ссылки. Ознакомлюсь.
  • +3
    Хайрер, Нёрсетт, Виннер, «Нежёсткие задачи» — СРОЧНО ЧИТАТЬ!

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

    Чем мне нравится один из этих методов Дормана-Принса — мы имеем каждый шаг в результате с пятым порядком точности (РК пятого порядка), но любые промежуточные точки можем рассчитать с четвёртым порядком, если сохранили промежуточные результаты вычислений.

    Я как-то развлекался, писал код на C для интегрирования систем из порядка 10000 уравнений, чтобы построить складки, которые образует, например, тонкий квадратныq резиновый коврик, опущеный на шар. Рендерил всё это в Pov-ray. В основном так и работал — ДП 5(4) с автокоррекцией.
  • 0
    Мне однажды в универе задали задачку — сделать модель солнечной системы на паскале. Я сделал, физику всю сохранил. Однако динамика была реальная и смотреть было на эту вселенную крайне скучно — она фактически стояла на месте.
    Когда я попытался внести в формулу сжатие времени для ускорения, луна оторвалась от земли, полетела в сторону марса, но промазала. Сделав круг вокруг солнца набрала 3ю космическую скорость и покинула солнечную систему. Такая вот грустная история, победить эту проблему так и не удалось.
    • +1
      Вы случайно не допустили классическую ошибку для подобных задач — прорисовка каждого полученного момента времени? Дело в том, что сама прорисовка значительно медленнее численного интегрирования, поэтому надо было не менять формулу, а просто рисовать каждый 100й кадр, пропуская остальные 99.
      • 0
        :-)
  • 0
    ЧЕРТ ПОБЕРИ!

    У меня получилось!
    Несколько дней изучения и чтения «Нежестких задач» и я осилил метод Рунге-Кутты!!!
    Мне удалось сбалансировать систему в новом проекте за 5 минут, при том, что методом Эйлера я целую неделю мучался подбором параметров и стабилизацией траектории!!!

    Всем огромное спасибо!
    • 0
      У меня в первой версии тоже вроде всё стабильно было со спутниками (правда, проверка была не очень длительная, и только в системе звезда–планет–спутник)

      Кстати, чтобы орбита была круговая, нужно чтобы скорость спутника была равна второй космической скорости по отношению с планете:
      speed = math.sqrt(2.0 * planet.mass / math.hypot(planet.x - satellite.x, planet.y - satellite.y))
      

      (это если G=1, как в данном симуляторе).
      Например, если в конфиге местоположение для Земли и Луны (соотношение масса 81:1) прописаны как
      [Earth]
      X=50
      Y=370
      Mass=81
      
      [Moon]
      X=30
      Y=370
      Mass=1
      

      и для Земли скорость
      VX=0.0
      VY=5.0
      

      то для Луны она будет
      VX = Earth.VX = 0
      VY = Earth.VY + sqrt(2.0 * Earth.mass / hypot(Earth.X - Moon.X, Earth.Y - Moon.Y)) = 8.16
      

      • 0
        Окститесь, вторая космическая по отношению к данному объекту — это выход на параболическую траекторию ухода на бесконечность от данного объекта.
        • 0
          Ох глюканул, сорри, первая конечно же имелась ввиду… Но тогда двойка лишняя в уравнении… странно, у меня именно в таком виде как записал получалась четко–круговая орбита. Проверю…
          • 0
            Проверил. Так и есть – двойки не надо. правильно будет
            speed = math.sqrt(planet.mass / math.hypot(planet.x - satellite.x, planet.y - satellite.y))
            
    • 0
      Терзают меня смутные сомнения…

      Я записал сумму всех сил, разложил по осям Х и У
      И для каждой из осей применил метод Рунге-Кутты.

      Соменение состоит в том, не потерял ли я чего-то важного при таком подходе. Все-же координаты зависимы, т.е. делать так нельзя. Но с другой стороны эта связь учтена в формуле:
              ax = (сумма) M0 * (X0 - x) / r**3
              ay = (сумма) M0 * (Y0 - y) / r**3
      

      где r — расстояние между двумя телами.

      Спасибо всем кто подскажет.
      • 0
        В данном случае, свойства системы таковы, что переменные смешиваются только в знаменателе под корнем. Так что, вроде бы правильно всё.

        Перечитайте habrahabr.ru/post/197754/#comment_6860970 — может быть, поймёте общую идею :)

        Я уже задумывался над постом о решении таких задач. Если удастся, до я доведу это до конца ;) (Ну, попробую написать то же самое, что в комментарии по ссылке, но понятным языком.)
        • 0
          :-)

          Спасибо за веру в мои скромные способности.

          Этот комментарий и послужил толчком для переделки :-)

          Ясно. Я провел некоторые исследования для системы из двух тел. Действительно, траектория остается на своем месте.
          Вот система из двух тел, расстояние между ними при наименьшем и наибольшем удалении:

          DIST MAX:  399.9574673146623
          DIST MIN:  76.50692564029397
          DIST MAX:  400.05701238567127
          DIST MIN:  76.29682103619999
          DIST MAX:  400.10700291341266
          DIST MIN:  76.61843705129094
          DIST MAX:  400.1784335186475
          DIST MIN:  76.31096590098528
          DIST MAX:  400.2307990384789
          DIST MIN:  76.5429266917722
          DIST MAX:  400.30332627343995
          DIST MIN:  76.22592172937536
          DIST MAX:  400.3598982590168
          DIST MIN:  76.33347805866529
          DIST MAX:  400.41598111793627
          DIST MIN:  76.44551804585252
          DIST MAX:  400.47259940381394
          DIST MIN:  76.53783120028531
          DIST MAX:  400.53049654382465
          DIST MIN:  76.5940516130298
          


          Еще раз спасибо.

Только зарегистрированные пользователи могут оставлять комментарии. Войдите, пожалуйста.