Хабр Курсы для всех
РЕКЛАМА
Практикум, Хекслет, SkyPro, авторские курсы — собрали всех и попросили скидки. Осталось выбрать!
Хабр доступен 24/7 благодаря поддержке друзей
Комментарии 69
Вот так вот и рождаются игры типа Kerbal Space Program, Free Orion, Solar (2).
Чего нельзя изобразить
books.google.ru/books?id=i7V2LXYdC-0C&pg=PA145
books.google.ru/books?id=i7V2LXYdC-0C&pg=PA145
То же такие мысли были. Я правда сразу за 3D взялся, время от времени кодю. Но у меня ещё ничего не летает. А вам удачи в проекте)
Вам мой положительный отзыв: развивайте проект! Вы мне напомнили давнюю мою мечту, о которой я уже давно забыл, а сейчас она возродилась, так что буду следить.
А можно уже новую версию опробовать?
А можно уже новую версию опробовать?
Спасибо за отзыв. Новая версия тут:
github.com/variostudio/spacesim
Исполняемый файл: main.py
Есть ключ -h
Есть ключ -f для загрузки .ini файла с конфигурацией системы
Есть несколько конфигураций. Пробуйте. Пишите что получилось.
Возможно, что-то криво закоммичено. Пишите, вечером поправлю
github.com/variostudio/spacesim
Исполняемый файл: main.py
Есть ключ -h
Есть ключ -f для загрузки .ini файла с конфигурацией системы
Есть несколько конфигураций. Пробуйте. Пишите что получилось.
Возможно, что-то криво закоммичено. Пишите, вечером поправлю
Не получилось поставить ConfigParser почему–то (Win7/Python2.6), но развил версию приведёную здесь до двух планет и звезда тоже обращается вокруг общего центра масс. Захватывает… Глядишь, разовьётся в вот такое:
Досадно…
Я, как начинающий в Питоне, писал сразу под версию 3.3
Может, в этом проблема?
Теоритически вы можете заменить ConfigParser на что-то из набора Python2.6
Там есть несколько готовых систем.
Я, как начинающий в Питоне, писал сразу под версию 3.3
Может, в этом проблема?
Теоритически вы можете заменить ConfigParser на что-то из набора Python2.6
Там есть несколько готовых систем.
Не очень понятно – тоже новичёк.
Кстати, чтобы поделится программой с теми кто не имеет Питона и сидит на виндах, можете скомпиливать в exe с помощью py2exe, а лучше с pywin32 (но может это только из–под винды будет работать...)
П.С. Планет уже 2, полная нестабильность, пока одну не вышвырнет из системы… или все разлетаются. Такие пируэты выдают! Как вы и писали, вся фишка в начальных данных, можно поймать очень красивую комбинацию.
Кстати, чтобы поделится программой с теми кто не имеет Питона и сидит на виндах, можете скомпиливать в exe с помощью py2exe, а лучше с pywin32 (но может это только из–под винды будет работать...)
П.С. Планет уже 2, полная нестабильность, пока одну не вышвырнет из системы… или все разлетаются. Такие пируэты выдают! Как вы и писали, вся фишка в начальных данных, можно поймать очень красивую комбинацию.
Я попробую откомпилировать вечером в.ехе
Да, иногда бывает сложно поймать стабильную комбинацию :-) Чем больше планет тоем тоньше грань.
Попробуйте Питон 3.3, если есть возможность.
Все должно запуститься.
Да, иногда бывает сложно поймать стабильную комбинацию :-) Чем больше планет тоем тоньше грань.
Попробуйте Питон 3.3, если есть возможность.
Все должно запуститься.
Ага. Прочитал, что до 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
Может, это поможет.
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
Может, это поможет.
Получилось с Питоном версии 3.2, версия 3.3 не вязалась с pygame. Спасибо!
Удалось разные конфигурации запустить (ключ -f и имя .ini файла)?
Да, всё отлично работает!
Кстати, в своей версии я выводил ещё текстовую информацию о расстояниях в окно. Возможно, кому-то будет интересно наблюдать как меняются параметры орбиты той планеты, на которую щёлкнули мышью. Также, интересно поэкспериментировать с тяжёлыми планетами, но тогда нужен режим (с возможностью отключить) где звезда не стационарна, а тоже динамична как и планеты. Я уже написал, что в своей версии я запускал планеты с массами близкими к звёздной — красота :)
Кстати, в своей версии я выводил ещё текстовую информацию о расстояниях в окно. Возможно, кому-то будет интересно наблюдать как меняются параметры орбиты той планеты, на которую щёлкнули мышью. Также, интересно поэкспериментировать с тяжёлыми планетами, но тогда нужен режим (с возможностью отключить) где звезда не стационарна, а тоже динамична как и планеты. Я уже написал, что в своей версии я запускал планеты с массами близкими к звёздной — красота :)
Видео есть? Если оно летает уже
Делал, что-то подобное на js + canvas. Планеты у меня можно добавлять кликом и потянуть для задания размера.
А как задавались начальные скорости для планет?
Рандомом, сейчас чуток переделал на ваши формулы, чуток изменю логику.
В начальных скоростях вся соль. Все великие столкновения планет или побеги из системы — суть начальных скоростей. Рандомом нельзя…
Насчет логики. Сначала нужно пересчитать ВСЕ ускорения а ТОЛЬКО ПОТОМ — все скорости и координаты.
Нельзя считать одновременно ускорения, скорости и координаты в цикле по объктам, иначе новые значения координат будут влиять на ускорения других объектов.
Об этом будет вторая статья.
Насчет логики. Сначала нужно пересчитать ВСЕ ускорения а ТОЛЬКО ПОТОМ — все скорости и координаты.
Нельзя считать одновременно ускорения, скорости и координаты в цикле по объктам, иначе новые значения координат будут влиять на ускорения других объектов.
Об этом будет вторая статья.
спасибо, а то и правда сумашествие какое-то происходит
Насчет логики. Сначала нужно пересчитать ВСЕ ускорения а ТОЛЬКО ПОТОМ — все скорости и координаты.
Нельзя считать одновременно ускорения, скорости и координаты в цикле по объктам, иначе новые значения координат будут влиять на ускорения других объектов.
Ну вы вот сами запутались.
Законы Ньютона для N тел напрямую дают систему из 3N дифференциальных уравнений второго порядка вида:
Здесь разные x — это какая-то из трёх координат каждого тела, всего 3N координат — итого 3N таких уравнений. x — координаты, x' — скорости, x'' — ускорения (ну, т.е. скорости и ускорения — первые и вторые производные координат)/
Функции в правой части на самом деле не зависят от скоростей и времени, но я просто хотел записать общий случай. Они, собственно говоря, сумма гравитационных взаимодействий, зависят только от координат.
Каждое такое уравнение эквивалентно двум уравнениям первого порядка, а вся система эквивалентна системе из 6N уравнений первого порядка вида (ещё я учёл, что гравитация от скорости и времени не зависит):
(т.е. 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 скоростей}, то вы храните предыдущие координаты и скорости. Хранить и пересчитывать ускорения вообще не нужно — они после дискретизации задачи (перехода к конечным разностям) вообще пропадают из всех формул. (Если строго говорить, их в формулах и не было — только в виде второй производной от координаты или первой от скорости.)
Да. Я именно про п.3 говорил.
Попервах у меня наслаивалось следующее значение на предыдущее. Простой баг был :-)
Попервах у меня наслаивалось следующее значение на предыдущее. Простой баг был :-)
Запишу задачу так:
X — это тот самый вектор
А это — уравнение, которое мы решаем. Для двух тел правая часть имеет вид (от t она явно не зависит):
%20%3D%20%5Cbegin%7Bbmatrix%7D%20k%20m_2%20%28x_1-x4%29%20%5C%5C%20k%20m_2%20%28x_2-x_5%29%20%5C%5C%20k%20m_2%20%28x_3-x_6%29%20%5C%5C%20k%20m_1%20%28x_4-x_1%29%20%5C%5C%20k%20m_1%20%28x_5-x_2%29%20%5C%5C%20k%20m_1%20%28x_6-x_3%29%20%5C%5C%20x_1%20%5C%5C%20..%20%5C%5C%20x_6%20%5Cend%7Bbmatrix%7D.gif)
где k — это общий знаменатель
Мог напутать разве что со знаками :) И забыл гравитационную постоянную, конечно же. Входит в k в виде ещё одного сомножителя.
x1..x3 — координаты первого тела, x4..x6 — второго тела, x7..x9 — скорости первого, x10..x12 — второго.
X — это тот самый вектор
А это — уравнение, которое мы решаем. Для двух тел правая часть имеет вид (от t она явно не зависит):
где k — это общий знаменатель
Мог напутать разве что со знаками :) И забыл гравитационную постоянную, конечно же. Входит в k в виде ещё одного сомножителя.
x1..x3 — координаты первого тела, x4..x6 — второго тела, x7..x9 — скорости первого, x10..x12 — второго.
Если вам интересно, то могу написать пост, как смоделировать закон Кулона (ну или всемирного тяготения — просто константы другие :) ). Я когда-то писал на Хабре один комментарий на эту тему, но народ потом переключился на обсуждение физики, решение диф. уравнений и так далее.
Посмотрите, там есть две старые видюшки моего JavaScript + Canvas решения. Скорей всего то, что d4rkr00t делал давно.
Посмотрите, там есть две старые видюшки моего JavaScript + Canvas решения. Скорей всего то, что d4rkr00t делал давно.
вообще да, было бы очень интересно, у меня то там совсем все просто было codepen.io/d4rkr00t/pen/hItgE движение по окружности + что-то типа силы грацитации.
Спасибо за комментарий и ссылку.
Да, очень интересно. Был бы рад увидеть Вашу статью.
И на программу d4rkr00t тоже интересно взглянуть :-)
Да, очень интересно. Был бы рад увидеть Вашу статью.
И на программу d4rkr00t тоже интересно взглянуть :-)
Лет 5 назад пользовался подобной программой, в ней тела при столкновении слипались, и можно было симулировать образование планет. Неплохо бы подобную вещь, но со сбалансированным слипанием, зависящим от притяжения.
Добрый день. Несколько замечаний по алгоритму.
Закон всемироного тяготения у вас записан не совсем верно: во-первых, это сила, во-вторых у нее две массы (массы обоих тел). Конечно, если силу подставить во второй закон Ньютона и расписать на оси, так и получится, но формулировка у вас дана неверная.
Отсюда вывод — сделайте его как силу, а в уравнения движения (которые в вашем случае являются интегрированием методом Эйлера, который имеет только одним преимуществом — простотой, но недостатков у него вагон) подставить силу, деленную на собственную массу. Это позволит улучшить модель и сделать взаимное влияние планет друг на друга (а оно очень большое), да и можно будет просто добавить, например, Луну или другие спутники. В идеале это еще позволит учитывать кучу возмущений, но вам вряд ли пока это нужно.
Второе — познакомьтесь с кеплеровскими элементами орбиты. Опять две причины: а — они позволяют не париться о выборе начальных условий (собственно, вы задаете элементы орбиты, а дальше они пересчитываются в XYZ, vx, vy, vz), б — в небесной механики координаты тел задаются как раз ими, так что если захотите брать реальные цифры — обязательно с ними столкнетесь.
Закон всемироного тяготения у вас записан не совсем верно: во-первых, это сила, во-вторых у нее две массы (массы обоих тел). Конечно, если силу подставить во второй закон Ньютона и расписать на оси, так и получится, но формулировка у вас дана неверная.
Отсюда вывод — сделайте его как силу, а в уравнения движения (которые в вашем случае являются интегрированием методом Эйлера, который имеет только одним преимуществом — простотой, но недостатков у него вагон) подставить силу, деленную на собственную массу. Это позволит улучшить модель и сделать взаимное влияние планет друг на друга (а оно очень большое), да и можно будет просто добавить, например, Луну или другие спутники. В идеале это еще позволит учитывать кучу возмущений, но вам вряд ли пока это нужно.
Второе — познакомьтесь с кеплеровскими элементами орбиты. Опять две причины: а — они позволяют не париться о выборе начальных условий (собственно, вы задаете элементы орбиты, а дальше они пересчитываются в XYZ, vx, vy, vz), б — в небесной механики координаты тел задаются как раз ими, так что если захотите брать реальные цифры — обязательно с ними столкнетесь.
Спасибо за комментарий.
Сейчас я пишу вторую часть статьи, где учтено взаимное влияние тел друг на друга. В примерах есть симуляция системы «Звезда — планета — спутник». Даже в двух вариантах.
Можно ли будет обратиться к Вам для для комментария обновленной версии?
Кстати, вопрос. В новой версии я пользуюсь этим же алгоритмом, но вычисляемом для каждой пары тел в обе стороны.
Т.е. есть тела 1...N
Я вычисляю ускорение при действии тела №i на тело №j, потом наоборот, тела №j на тело №i
Все суммируется. Получается N^2 вычислений для системы из N тел. Каждое взаимодействие учтено.
Этот подход имеет право на существование?
Еще раз спасибо за конструктивный комментарий.
Сейчас я пишу вторую часть статьи, где учтено взаимное влияние тел друг на друга. В примерах есть симуляция системы «Звезда — планета — спутник». Даже в двух вариантах.
Можно ли будет обратиться к Вам для для комментария обновленной версии?
Кстати, вопрос. В новой версии я пользуюсь этим же алгоритмом, но вычисляемом для каждой пары тел в обе стороны.
Т.е. есть тела 1...N
Я вычисляю ускорение при действии тела №i на тело №j, потом наоборот, тела №j на тело №i
Все суммируется. Получается N^2 вычислений для системы из N тел. Каждое взаимодействие учтено.
Этот подход имеет право на существование?
Еще раз спасибо за конструктивный комментарий.
Не очень понял с ускорениями.
Я бы советовал двигаться в сторону второго закона Ньютона. А именно — ускорение считать как 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, но лучше засекать время с предыдущей итерации).
Да, ко мне можно обращаться, всегда с радостью помогу.
Я бы советовал двигаться в сторону второго закона Ньютона. А именно — ускорение считать как 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, но лучше засекать время с предыдущей итерации).
Да, ко мне можно обращаться, всегда с радостью помогу.
Не очень понял с ускорениями.
Я бы советовал двигаться в сторону второго закона Ньютона. А именно — ускорение считать как Ai = Sum(Fij)/m, где Fij — сила взаимного притяжения тел (Fij = G * mi * mj * (Xi — Xj) / (rij ^ 3), тут записано в векторной форме).
Я ж точно так и делал. Только сразу сократил m и оперировал не понятием «сила», а понятием «ускорение» :-)
Мне уже подсказали воспользоваться другими численными методами. Возможно, я так и поступлю, если удастся осилить. Университетский курс численных методов требует много времени. Школьный курс физики намного менее требовательный.
Я бы советовал двигаться в сторону второго закона Ньютона. А именно — ускорение считать как Ai = Sum(Fij)/m, где Fij — сила взаимного притяжения тел (Fij = G * mi * mj * (Xi — Xj) / (rij ^ 3), тут записано в векторной форме).
Я ж точно так и делал. Только сразу сократил m и оперировал не понятием «сила», а понятием «ускорение» :-)
Мне уже подсказали воспользоваться другими численными методами. Возможно, я так и поступлю, если удастся осилить. Университетский курс численных методов требует много времени. Школьный курс физики намного менее требовательный.
Вообще хорошо бы вам почитать первый том Ландау-Лифшица. Сложно, да. Но там научат строить строго систему уравнений движения, дифференциальных, в данном случае они будут второго порядка. Ну да из системы N второго порядка сделать 2N первого — раз плюнуть. Это будет универсальный метод для любых количеств любых тел.
А эти уравнения можно интегрировать прямо сразу численно.
А эти уравнения можно интегрировать прямо сразу численно.
Использованный метод численного дифференцирования фактически представляет собой метод Эйлера первого порядка, и, в случае эллиптических орбит, будет быстро набегать систематическая ошибка — планеты будут потихоньку съезжать «от Солнца».
Традиционно в программах такого рода используют метод Рунге-Кутты 4-го порядка. Он лишь чуть сложнее метода Эйлера в программировании, но намного точнее.
Традиционно в программах такого рода используют метод Рунге-Кутты 4-го порядка. Он лишь чуть сложнее метода Эйлера в программировании, но намного точнее.
Точно, так и есть…
Буду осваивать новый метод. Есть над чем поработать теперь.
Огромное спасибо.
Буду осваивать новый метод. Есть над чем поработать теперь.
Огромное спасибо.
Дополню. Во-первых, используется не сам метод Ронге-Кутта, а его модификация — Дормана-Принса. Он хорош тем, что пересчитывает шаг и где можно шагать без потерь огромными шагами так и делает, а где нужно точно проинтегрировать — уменьшает шаг. Во-вторых, очень удобно в подобных, да и во всех вообще случаях. Но если вам вдруг взбредет симулировать движение планет на продолжительный период времени, стоит обратить внимание на орбито-устойчивые методы. Они позволяют интегрировать сотни, тысячи, миллионы оборотов с достаточной точностью. С ходу ссылку найти не удалось, но могу подсказать литературу, если нужно.
ДП — няшечка. Обычно гоняют 5(4), но есть он и более высокого порядка — 8(7)
У меня никогда не хватало терпения вбить всю таблицу для 8(7) порядка))
Я вот давно хочу еще попробовать метод Адамса, но руки не доходят. Сейчас буду комплексировать GPS/ГЛОНАСС и радио-баровысотомер (модели, конечно), попробую написать Адамса или еще что-то из экзотики)
Я вот давно хочу еще попробовать метод Адамса, но руки не доходят. Сейчас буду комплексировать GPS/ГЛОНАСС и радио-баровысотомер (модели, конечно), попробую написать Адамса или еще что-то из экзотики)
+1 про ДП8(7) :) я по-моему его вбил, считало неправильно, искать ошибку… хрен найдёшь её.
Есть ещё неявные методы семейства РК, там свои плюсы (устойчивее) и минусы (реализовывать посложней будет).
Адамс хорош, когда реально должно получаться что-то очень близкое к многочленам. Он же аппроксимирует многочленами. Если точно решение — многочлен порядка K, то метод Адамса любоого порядка P>K должен давать 0 в старших членах, т.к. он точно найдёт K+1 коэффициентов этого многочлена. Отсюда пляшем выбор метода интегрирования.
Есть ещё неявные методы семейства РК, там свои плюсы (устойчивее) и минусы (реализовывать посложней будет).
Адамс хорош, когда реально должно получаться что-то очень близкое к многочленам. Он же аппроксимирует многочленами. Если точно решение — многочлен порядка K, то метод Адамса любоого порядка P>K должен давать 0 в старших членах, т.к. он точно найдёт K+1 коэффициентов этого многочлена. Отсюда пляшем выбор метода интегрирования.
За РК четвёртого — минус вам. Ну хватит уже его насиловать. Морально устарел РК 4-го порядка, никто его не использует, только некоторые, кто «что-то слышал».
Когда нам удастся разобрать этот случай, мы легко перейдем к сложным системам со взаимным влиянием звезд и планет друг на друга.
На этой фразе наверное все физики рассмеялись.
В рамках методов прямого численного дифференцирования N-body систем, как у автора — действительно легко.
Вот если аналитически — тут да, уже решение задачи трёх тел весьма проблематично.
Вот если аналитически — тут да, уже решение задачи трёх тел весьма проблематично.
Да нууууу? Легко, да?
Как вы думаете, в системе из 10 тел, далеко вы уйдёте методом Эйлера, чтоб погрешности не выросли хотя бы на 10%?
Очень недалеко.
Как вы думаете, в системе из 10 тел, далеко вы уйдёте методом Эйлера, чтоб погрешности не выросли хотя бы на 10%?
Очень недалеко.
> Чтобы решить задачу, нужно сначала четко себе ее представить.
Вот-вот.
Рекомендую ознакомится с исходниками давно построенного и проверенного велосипеда под названием "Симулятор Солнечной системы".
А именно, что б долго не искать, с классом Orbit и, например, его использованием в методе Body::getPosition().
p.s. Ваша заметка плохая. Теории нет, решения заявленной задачи нет. То, что вы накодили — кривое и, строго говоря, нерабочее… Вообще оставляет впечатление «нам вчера в классе рассказали про питон и зацените, что я сегодня набросал».
Не надо такую же вторую писать.
Вот-вот.
Рекомендую ознакомится с исходниками давно построенного и проверенного велосипеда под названием "Симулятор Солнечной системы".
А именно, что б долго не искать, с классом Orbit и, например, его использованием в методе Body::getPosition().
p.s. Ваша заметка плохая. Теории нет, решения заявленной задачи нет. То, что вы накодили — кривое и, строго говоря, нерабочее… Вообще оставляет впечатление «нам вчера в классе рассказали про питон и зацените, что я сегодня набросал».
Не надо такую же вторую писать.
Хайрер, Нёрсетт, Виннер, «Нежёсткие задачи» — СРОЧНО ЧИТАТЬ!
Там вы найдёте не только коэффициенты методов Рунге-Кутты, но и узнаете про методы Адамса (уход несколько в другую сторону, нежели коррекция за счёт всё больших степеней производных), а также методы а-ля Дормана-Принса с автоматической коррекцией величины шага и оценкой погрешности при интегрировании.
Чем мне нравится один из этих методов Дормана-Принса — мы имеем каждый шаг в результате с пятым порядком точности (РК пятого порядка), но любые промежуточные точки можем рассчитать с четвёртым порядком, если сохранили промежуточные результаты вычислений.
Я как-то развлекался, писал код на C для интегрирования систем из порядка 10000 уравнений, чтобы построить складки, которые образует, например, тонкий квадратныq резиновый коврик, опущеный на шар. Рендерил всё это в Pov-ray. В основном так и работал — ДП 5(4) с автокоррекцией.
Там вы найдёте не только коэффициенты методов Рунге-Кутты, но и узнаете про методы Адамса (уход несколько в другую сторону, нежели коррекция за счёт всё больших степеней производных), а также методы а-ля Дормана-Принса с автоматической коррекцией величины шага и оценкой погрешности при интегрировании.
Чем мне нравится один из этих методов Дормана-Принса — мы имеем каждый шаг в результате с пятым порядком точности (РК пятого порядка), но любые промежуточные точки можем рассчитать с четвёртым порядком, если сохранили промежуточные результаты вычислений.
Я как-то развлекался, писал код на C для интегрирования систем из порядка 10000 уравнений, чтобы построить складки, которые образует, например, тонкий квадратныq резиновый коврик, опущеный на шар. Рендерил всё это в Pov-ray. В основном так и работал — ДП 5(4) с автокоррекцией.
Мне однажды в универе задали задачку — сделать модель солнечной системы на паскале. Я сделал, физику всю сохранил. Однако динамика была реальная и смотреть было на эту вселенную крайне скучно — она фактически стояла на месте.
Когда я попытался внести в формулу сжатие времени для ускорения, луна оторвалась от земли, полетела в сторону марса, но промазала. Сделав круг вокруг солнца набрала 3ю космическую скорость и покинула солнечную систему. Такая вот грустная история, победить эту проблему так и не удалось.
Когда я попытался внести в формулу сжатие времени для ускорения, луна оторвалась от земли, полетела в сторону марса, но промазала. Сделав круг вокруг солнца набрала 3ю космическую скорость и покинула солнечную систему. Такая вот грустная история, победить эту проблему так и не удалось.
ЧЕРТ ПОБЕРИ!
У меня получилось!
Несколько дней изучения и чтения «Нежестких задач» и я осилил метод Рунге-Кутты!!!
Мне удалось сбалансировать систему в новом проекте за 5 минут, при том, что методом Эйлера я целую неделю мучался подбором параметров и стабилизацией траектории!!!
Всем огромное спасибо!
У меня получилось!
Несколько дней изучения и чтения «Нежестких задач» и я осилил метод Рунге-Кутты!!!
Мне удалось сбалансировать систему в новом проекте за 5 минут, при том, что методом Эйлера я целую неделю мучался подбором параметров и стабилизацией траектории!!!
Всем огромное спасибо!
У меня в первой версии тоже вроде всё стабильно было со спутниками (правда, проверка была не очень длительная, и только в системе звезда–планет–спутник)
Кстати, чтобы орбита была круговая, нужно чтобы скорость спутника была равна второй космической скорости по отношению с планете:
(это если G=1, как в данном симуляторе).
Например, если в конфиге местоположение для Земли и Луны (соотношение масса 81:1) прописаны как
и для Земли скорость
то для Луны она будет
Кстати, чтобы орбита была круговая, нужно чтобы скорость спутника была равна второй космической скорости по отношению с планете:
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
Окститесь, вторая космическая по отношению к данному объекту — это выход на параболическую траекторию ухода на бесконечность от данного объекта.
Терзают меня смутные сомнения…
Я записал сумму всех сил, разложил по осям Х и У
И для каждой из осей применил метод Рунге-Кутты.
Соменение состоит в том, не потерял ли я чего-то важного при таком подходе. Все-же координаты зависимы, т.е. делать так нельзя. Но с другой стороны эта связь учтена в формуле:
где r — расстояние между двумя телами.
Спасибо всем кто подскажет.
Я записал сумму всех сил, разложил по осям Х и У
И для каждой из осей применил метод Рунге-Кутты.
Соменение состоит в том, не потерял ли я чего-то важного при таком подходе. Все-же координаты зависимы, т.е. делать так нельзя. Но с другой стороны эта связь учтена в формуле:
ax = (сумма) M0 * (X0 - x) / r**3
ay = (сумма) M0 * (Y0 - y) / r**3
где r — расстояние между двумя телами.
Спасибо всем кто подскажет.
В данном случае, свойства системы таковы, что переменные смешиваются только в знаменателе под корнем. Так что, вроде бы правильно всё.
Перечитайте habrahabr.ru/post/197754/#comment_6860970 — может быть, поймёте общую идею :)
Я уже задумывался над постом о решении таких задач. Если удастся, до я доведу это до конца ;) (Ну, попробую написать то же самое, что в комментарии по ссылке, но понятным языком.)
Перечитайте habrahabr.ru/post/197754/#comment_6860970 — может быть, поймёте общую идею :)
Я уже задумывался над постом о решении таких задач. Если удастся, до я доведу это до конца ;) (Ну, попробую написать то же самое, что в комментарии по ссылке, но понятным языком.)
:-)
Спасибо за веру в мои скромные способности.
Этот комментарий и послужил толчком для переделки :-)
Ясно. Я провел некоторые исследования для системы из двух тел. Действительно, траектория остается на своем месте.
Вот система из двух тел, расстояние между ними при наименьшем и наибольшем удалении:
Еще раз спасибо.
Спасибо за веру в мои скромные способности.
Этот комментарий и послужил толчком для переделки :-)
Ясно. Я провел некоторые исследования для системы из двух тел. Действительно, траектория остается на своем месте.
Вот система из двух тел, расстояние между ними при наименьшем и наибольшем удалении:
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
Еще раз спасибо.
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Создаем симулятор солнечной системы