Создаем симулятор солнечной системы

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

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

спасибо, а то и правда сумашествие какое-то происходит

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

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

Законы Ньютона для 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 скоростей}, то вы храните предыдущие координаты и скорости. Хранить и пересчитывать ускорения вообще не нужно — они после дискретизации задачи (перехода к конечным разностям) вообще пропадают из всех формул. (Если строго говорить, их в формулах и не было — только в виде второй производной от координаты или первой от скорости.)

старая версию по ссылке выше, на codepen.io

У меня никогда не хватало терпения вбить всю таблицу для 8(7) порядка))
Я вот давно хочу еще попробовать метод Адамса, но руки не доходят. Сейчас буду комплексировать GPS/ГЛОНАСС и радио-баровысотомер (модели, конечно), попробую написать Адамса или еще что-то из экзотики)

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

И это — всего одна переменная.

Проверил. Так и есть – двойки не надо. правильно будет

speed = math.sqrt(planet.mass / math.hypot(planet.x - satellite.x, planet.y - satellite.y))