Телескопы покупают здесь


A A A A Автор Тема: Упрощенное моделирование движения в гравитационном поле  (Прочитано 1667 раз)

0 Пользователей и 1 Гость просматривают эту тему.

Оффлайн Небесный МеханикАвтор темы

  • Новичок
  • *
  • Сообщений: 20
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Небесный Механик
Пишу программу для моделирования движения межпланетного космического корабля. Двигатель корабля может быть включен сколь угодно долго, поэтому никакого уравнения траектории нет, и имеет смысл пошаговое моделирование.

Опишу алгоритм, используемый на данный момент:
1. Имеются начальные координаты (x,y,z) и вектор скорости по осям: v(x,y,z)
2. Вычисляется смещение за одну секунду, как если бы оно происходило с постоянной скоростью v. Получаем координаты x1,y1,z1
3. В новой точке вычисляется векторная сумма сил, действующих на тело (от гравирующих объектов и от двигателей). Затем получаем вектор ускорения: a(x1,y1,z1)
4. Вычисляем новый вектор скорости: v(x1,y1,z1) = v(x,y,z) + a(x1,y1,z1)
5. Повторяем итерацию.

Проверку данного алгоритма провел на нескольких планетах, чтобы оценить точность метода. Задавая в качестве начальных данных расстояние и скорость в перигелии. Результаты получились в целом неплохие:
1. Смоделированные орбиты визуально неотличимы от эллиптических.
2. Расстояние и скорость в афелии определяется с высокой точностью.
3. Период на доли процента отличается от реального: для Земли 367 дней, для Марса также около 2 дней дольше.

Баг с периодом, в принципе, не критичен, но хотелось бы узнать, существует ли простой способ повысить точность.

Вопросы такие:
- Правильно ли я понимаю, что по сути провожу численное интегрирование методом прямоугольника?
- Существует ли простой способ перейти от метода прямоугольника к методу трапеций, то есть на каждом шаге делать аппроксимацию ускорения и скорости неким серединным значением? И повысит ли это точность периода?

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

Оффлайн xd

  • *****
  • Сообщений: 17 977
  • Благодарностей: 378
    • Skype - deimos.belastro.net
  • Награды Открытие комет, астероидов, сверхновых звезд, научно значимые исследования.
    • Сообщения от xd
    • Белорусская любительская астрономическая сеть
Предиктивные схемы численного интегрирования в помощь.
У природы нет плохой погоды, у неё просто на нас аллергия.

Учение без размышления бесполезно, но и размышление без учения опасно /Конфуций/
Слово есть поступок. /Л. Толстой/

Оффлайн Geen

  • *****
  • Сообщений: 12 210
  • Благодарностей: 200
  • Мне нравится этот форум!
    • Сообщения от Geen
Смотрите методы Рунге-Кутта http://en.wikipedia.org/wiki/Runge-Kutta_methods, адаптивные методы Рунге-Кутта, симплектические методы http://en.wikipedia.org/wiki/Symplectic_integrator (хороши тем, что время точно обратимо).
Если не пугает JavaScript, то можно посмотреть на моей странице http://www.k-labs.ru/dms/nbody.html примеры реализаций и ссылки на другие методы.
Могу отметить, что в "аналогичной программе" Orbiter http://orbit.medphys.ucl.ac.uk/ используется симплектический интегратор 6-го порядка по умолчанию.
Если у тебя есть фонтан, заткни его, дай отдохнуть и фонтану.

А ещё мы любим обсуждать вкус устриц с теми кто их ел...

Оффлайн Aluminium

  • *****
  • Сообщений: 2 619
  • Благодарностей: 99
  • Мне нравится этот форум!
    • Сообщения от Aluminium
3. Период на доли процента отличается от реального: для Земли 367 дней, для Марса также около 2 дней дольше.

Баг с периодом, в принципе, не критичен, но хотелось бы узнать, существует ли простой способ повысить точность.
Дело может быть в неточных начальных данных. Или гравитационных параметрах тел.

Способы повысить точность попроще - это либо модифицированный метод Эйлера, либо чуть посложнее метод Рунге-Кутты 4-го порядка.

Оффлайн Небесный МеханикАвтор темы

  • Новичок
  • *
  • Сообщений: 20
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Небесный Механик
Если не пугает JavaScript, то можно посмотреть на моей странице http://www.k-labs.ru/dms/nbody.html примеры реализаций и ссылки на другие методы.
Не удалось, говорят, получить доступ к сайту.

Оффлайн Небесный МеханикАвтор темы

  • Новичок
  • *
  • Сообщений: 20
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Небесный Механик
Смотрите методы Рунге-Кутта http://en.wikipedia.org/wiki/Runge-Kutta_methods, адаптивные методы Рунге-Кутта, симплектические методы http://en.wikipedia.org/wiki/Symplectic_integrator (хороши тем, что время точно обратимо).
Способы повысить точность попроще - это либо модифицированный метод Эйлера, либо чуть посложнее метод Рунге-Кутты 4-го порядка.
Я если честно просто затрудняюсь на стадии получения исходного уравнения. Исходная формула F =GmM/R^2, а интегрируя дважды, нам требуется получить, по сути, R, то есть рекурсия выходит. Вот на этом у меня и затык  :(

Оффлайн Geen

  • *****
  • Сообщений: 12 210
  • Благодарностей: 200
  • Мне нравится этот форум!
    • Сообщения от Geen

Если не пугает JavaScript, то можно посмотреть на моей странице http://www.k-labs.ru/dms/nbody.html примеры реализаций и ссылки на другие методы.
Не удалось, говорят, получить доступ к сайту.
Вроде бы поправил.
Если у тебя есть фонтан, заткни его, дай отдохнуть и фонтану.

А ещё мы любим обсуждать вкус устриц с теми кто их ел...

Оффлайн Aluminium

  • *****
  • Сообщений: 2 619
  • Благодарностей: 99
  • Мне нравится этот форум!
    • Сообщения от Aluminium
Я если честно просто затрудняюсь на стадии получения исходного уравнения.
Математическая формулировка задачи N тел.
В итоге будет система уравнений вида dz/dt = f(z, t). Где z - вектор из 6N переменных, а именно координат и скоростей тел, t - время.

Оффлайн Toth

  • *****
  • Сообщений: 2 581
  • Благодарностей: 174
    • Сообщения от Toth
Занимался я этим. И Рунге-Кутты , и Эрмита, и Эверхарта пробовал.
Лучше всего получилось по методу Эрмита ( там, где вычисляется рывок - производная ускорения ). Эверхарта - не получилось, у него точность у меня не выше, чем у Рунге-Кутта . Наверное, у меня руки кривые, ибо Эверхарта считается одним из крутейших.

Вот отсюда брал формулы - https://www.google.ru/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0ahUKEwjazMXJ_aTXAhWMIpoKHdASB-QQFggnMAA&url=http%3A%2F%2Fvadimchazov.narod.ru%2Ftext_pdf%2Freasint.doc&usg=AOvVaw1R4hWZYSXDlWXK8yw7yaCe
Только там опечатки есть. Вот - исправления

Оффлайн Небесный МеханикАвтор темы

  • Новичок
  • *
  • Сообщений: 20
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Небесный Механик
Вот отсюда брал формулы - https://www.google.ru/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0ahUKEwjazMXJ_aTXAhWMIpoKHdASB-QQFggnMAA&url=http%3A%2F%2Fvadimchazov.narod.ru%2Ftext_pdf%2Freasint.doc&usg=AOvVaw1R4hWZYSXDlWXK8yw7yaCe
Спасибо. Правильно я понимаю, что метод Эйлера на 4-й странице - это примерно то, что описано у меня в стартовом посте?

Оффлайн Geen

  • *****
  • Сообщений: 12 210
  • Благодарностей: 200
  • Мне нравится этот форум!
    • Сообщения от Geen
Вот отсюда брал формулы - https://www.google.ru/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0ahUKEwjazMXJ_aTXAhWMIpoKHdASB-QQFggnMAA&url=http%3A%2F%2Fvadimchazov.narod.ru%2Ftext_pdf%2Freasint.doc&usg=AOvVaw1R4hWZYSXDlWXK8yw7yaCe
Спасибо. Правильно я понимаю, что метод Эйлера на 4-й странице - это примерно то, что описано у меня в стартовом посте?
Да.
Если у тебя есть фонтан, заткни его, дай отдохнуть и фонтану.

А ещё мы любим обсуждать вкус устриц с теми кто их ел...

Оффлайн ulitkanasklone

  • *****
  • Сообщений: 8 645
  • Благодарностей: 280
  • тихо,тихо ползи улитка по склону Фудзи..
    • Сообщения от ulitkanasklone
    • СПРАВОЧНАЯ ИНФОРМАЦИЯ ПО ОТО, СТАТЬИ И НЕКОТОРЫЕ ЗАДАЧИ

Если не пугает JavaScript, то можно посмотреть на моей странице http://www.k-labs.ru/dms/nbody.html примеры реализаций и ссылки на другие методы.
Не удалось, говорят, получить доступ к сайту.
Вроде бы поправил.
Если честно, я вашей таблице ничего не понял. Вы бы где-нибудь на стартовой странице сделали бы пояснение, что моделируете , какие параметры вводятся и что в итоге мы получим.
Мне нравится этот форум
моя страница:
http://антониум.рф/ОТО/GT.html

Оффлайн Небесный МеханикАвтор темы

  • Новичок
  • *
  • Сообщений: 20
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Небесный Механик
Сделал по Рунге-Кутта, протестил на Марсе. Отметил странную особенность: в зависимости от выбранного шага период скачет, например для 500 секунд - он получился 687 суток при 502 секундах - 685, при 503 - снова 687. И так далее циклически. В целом "гуляет" от 683 до 688.
Это нормально, или искать баг?

Оффлайн xd

  • *****
  • Сообщений: 17 977
  • Благодарностей: 378
    • Skype - deimos.belastro.net
  • Награды Открытие комет, астероидов, сверхновых звезд, научно значимые исследования.
    • Сообщения от xd
    • Белорусская любительская астрономическая сеть
В формулах ошибок метода. Там остаточная погрешность сидит.
У природы нет плохой погоды, у неё просто на нас аллергия.

Учение без размышления бесполезно, но и размышление без учения опасно /Конфуций/
Слово есть поступок. /Л. Толстой/

Оффлайн Aluminium

  • *****
  • Сообщений: 2 619
  • Благодарностей: 99
  • Мне нравится этот форум!
    • Сообщения от Aluminium
Шаг в 500 секунд мелковат, по-моему, для задачи про Марс и выбранного метода. Можно на порядок-другой увеличить. Вряд ли за 100 тысяч шагов может накопиться ошибка округления, если только Вы не используете одинарную точность. Проблема, скорее всего, в реализации метода или способе вычисления периода.

Оффлайн Toth

  • *****
  • Сообщений: 2 581
  • Благодарностей: 174
    • Сообщения от Toth
Я периоды не проверял. Я проверял саму траекторию.
Но у меня тоже получались колебания, то есть при меньшем шаге - чуть большая погрешность.
Вот - 2 примера.
1) Задача 2-х тел. ИСЗ на низкой почти круговой орбите (e=0.01) за 3 часа. Сравнивал с расчетами по Кеплеру.
2) Задача 10-ти тел - Солнечная система за 361 день . Тут - только по мет. Эрмита считал, РК для этого слишком тормозной. Исходные данные и контрольные - из http://vo.imcce.fr/webservices/miriade/?forms , модель DE\LE405
Обратите внимание на погрешность у Меркурия.

PS Попробую тоже период Марса посчитать в задаче просто 2-х тел.

Оффлайн Небесный МеханикАвтор темы

  • Новичок
  • *
  • Сообщений: 20
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Небесный Механик
Шаг в 500 секунд мелковат, по-моему, для задачи про Марс и выбранного метода. Можно на порядок-другой увеличить. Вряд ли за 100 тысяч шагов может накопиться ошибка округления, если только Вы не используете одинарную точность. Проблема, скорее всего, в реализации метода или способе вычисления периода.
Период я беру просто как разность во времени между прохождениями перигелия.
Ошибку я допускаю, но почему она не сказывается на значении афелия? 249232639 км, то есть очень близко к табличному. И в перигелии после первого оборота отклонение всего на 25 метров (не километров!) от значения, заданного в начальных условиях.

А период при этом гуляет. Не понимаю, короче, в чем дело.

Оффлайн Aluminium

  • *****
  • Сообщений: 2 619
  • Благодарностей: 99
  • Мне нравится этот форум!
    • Сообщения от Aluminium
Период я беру просто как разность во времени между прохождениями перигелия.
А как определяете прохождение перигелия? В задаче кроме Марса и Солнца больше ничего нет?

Может, проще посчитать положение на конкретный момент с разным шагом и посмотреть, как метод сходится? Только не 500 и 502 сравнивать, а различающиеся в два раза. Скажем, 25000, 50000, 100000, 200000. От уменьшения шага в два раза точность должна улучшаться на порядок. С совсем малым шагом начнут доминировать ошибки округления.

А период при этом гуляет. Не понимаю, короче, в чем дело.
Загадка природы.

Оффлайн Toth

  • *****
  • Сообщений: 2 581
  • Благодарностей: 174
    • Сообщения от Toth
Можете повторить по моим данным, и сравнить результаты.
Задача 2-х тел, Солнце и Марс. Данные координатам Марса взял условно, все равно - без остальных планет, просто идеальная невозмущенная орбита при a=1.523662 а.е.  и e = 0.0933941
Единицы у меня в системе СИ. Понимаю, что в небесной механике это - моветон, но можно перевести. ( а.е. = 149597870691 )
GM Солнца + GM  Марса = 1.3271244001799E+20 + 4.2828314258067E+13
Нач. вектор положения Марса =(206648658093.026;0;0), скорость =(0;26498.9066619695;0) - в перигелии. У Солнца - все 0.
Через 59353325.08 секунд Марс должен вернуться на прежнее место. В гелиоцентрических координатах ( то есть радиус-вектор Марса минус радиус-вектор Солнца. Хотя Марс легкий, но все же .. )
Значения шагов - не совсем точно те, потому что прога их подгоняет, чтобы было целое число шагов, потому что dT - не круглое число:
n:=Round(dT/h);
 if n=0 then n:=1;
 h:=dT/n;
Вот, результат - ошибка в конечном положении. У Эрмита - уже насыщение, там уже машинное эпсилон наверное сказывается.
« Последнее редактирование: 05 Ноя 2017 [03:07:27] от Toth »


Оффлайн Aluminium

  • *****
  • Сообщений: 2 619
  • Благодарностей: 99
  • Мне нравится этот форум!
    • Сообщения от Aluminium
Единицы у меня в системе СИ. Понимаю, что в небесной механике это - моветон, но можно перевести.
В принципе без разницы, в чём считать, лишь бы единицы были указаны. В системе СИ главное - гравитационную постоянную на массу не умножать.