Различия между числовыми пропагаторами

Я стажер, который работает над числовым орбитальным пропагатором, разработанным в компании. Я не могу показать вам код, но могу сказать, что распространитель был разработан для работы в Simulink. Моя работа заключалась в том, чтобы извлечь распространитель, чтобы использовать его в качестве скрипта Matlab. Чтобы оценить работу моего скрипта, я сравниваю полученные результаты с результатами другого числового распространителя (GMAT). Начальные элементы орбиты взяты с орбиты SSO Repeat Ground Track. В настоящее время я рассматриваю помехи, вызванные только J2, и я использую интегратор типа Рунге-Кутта четвертого порядка. Анализы с двумя пропагаторами за период 52 дня с фиксированным размером шага 60 секунд имеют разные тренды, показанные на рисунке (слева показаны элементы орбиты, а справа показан вектор состояния).

введите описание изображения здесь

Я не могу объяснить эту разницу, особенно в векторе состояния, из которого получаются орбитальные элементы. Сначала я подумал, что это связано с эталонными системами, поскольку данные GMAT находятся в ICRF, а данные моего пропагатора — в GCRS. Поэтому я использовал те же функции преобразования для данных GMAT, и при сравнении (этих новых результатов) с исходными не было ошибок, кроме тех, что связаны с числовыми. Так что я подумал, что ошиблась функция геопотенциала. Но даже в этом случае при прохождении одного и того же вектора состояния из GMAT ускорение, рассчитанное функцией, было таким же, как и предоставленное программой. Поэтому я думаю, что это может быть ошибка интегратора, так как есть периодический дрейф ошибки, показанный на следующем рисунке.введите описание изображения здесь

Однако Рунге-Кутта, по сравнению с тем, что я нашел в Интернете, кажется, настроен хорошо. Можете ли вы назвать причину, объясняющую эти расхождения?

Я добавляю больше графиков, чтобы показать вам тренд орбитальных элементов и вектора состояния за 24 часа. Со временем различия имеют тенденцию увеличиваться. В частности, мой пропагатор кажется ускоренным по сравнению с GMAT.

введите описание изображения здесь введите описание изображения здесь введите описание изображения здесь

Интересный вопрос и прямо по теме здесь, в Space SE!
Крайне важно знать: эти справочные данные GMAT рассчитаны с точно такой же установкой, гравитационным потенциалом и т. д., или просто с установкой, которая, по вашему мнению , должна быть такой же?
Да, я установил одинаковую степень и порядок для геопотенциала, который является единственным возмущением, которое я рассматриваю, а также те же элементы орбиты, эпоху и интегратор (RK4 с фиксированным размером шага 60 секунд). Единственная разница заключается в том, что GMAT рассматривает модель EGM-96, в то время как мой распространитель рассматривает EGM-2008 и основан на соглашениях IERS Tech Note 36. Также результаты GMAT находятся в ICRF, а мои — в GCRS. Именно по этой причине я ожидал других результатов, но не такой высокой ошибки.
...ну это уже много чего может пойти не так. Расхождения почти наверняка связаны с этими деталями, и решатель ОДУ не имеет к ним никакого отношения. Различные индексы/комплексы и т. д. для коэффициентов геопотенциала? Что-то не так с преобразованиями координат? И вы используете только квадруполь, а GMAT использует все?

Ответы (3)

Вы должны проверить, какой распространитель использует GMAT. Метод RK4 с постоянным размером шага со временем будет отклоняться от истинного решения ОДУ, поэтому, если GMAT использует метод с адаптивным размером шага, который в некоторых точках изменяется от порядка миллисекунд до десятков секунд, скорее всего, это будет ближе к решению, чем RK4.

По сути, это то, что делает постоянный размер шага RK4: введите описание изображения здесьон выполняет несколько оценок производной между настоящим и следующим временным шагом, чтобы получить лучшую оценку того, как функция изменяется (в отличие, скажем, от метода Эйлера 1-го порядка). ). Затем он берет средневзвешенное значение этих оценок, чтобы определить, какую производную он будет использовать для выполнения этого шага. Как видно из графика, в некоторых точках он занижает площадь под кривой, а в других завышает, что со временем приведет к возникновению ошибки.

Вот уменьшенная версия (и временной шаг 500 секунд), чтобы показать эффект.введите описание изображения здесь

Поскольку вы используете 60 секунд, ошибка будет расти медленнее, но, тем не менее, со временем будет дрейфовать.

Адаптивный решатель размера шага изменит размер шага в зависимости от того, насколько быстро меняется производная (вторая производная):введите описание изображения здесь

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

Наиболее распространенным из этих методов является RK4-5, который сравнивает решение методов RK4 и RK5 и в зависимости от того, насколько они различаются, соответствующим образом меняет временной шаг.

Таким образом, в целом с вашей моделью / реализацией RK4 может быть что-то не так, это может зависеть от того, какой распространитель GMAT использует, с которым вы сравниваете свое решение.

Также, если вам интересно, эти скриншоты взяты из видео, которое я сделал на решателях ODE:

@ChrisR Я не задавал вопрос? Я также не использую Simulink, Matlab или GMAT.
Привет @AlfonsoGonzalez, спасибо за ответ. В GMAT (чтобы иметь стабильные результаты) я использовал Runge Kutta 4 и заставил размер шага оставаться постоянным, установив 60 секунд в качестве минимального и максимального размера шага. По этой причине я ожидал получить те же результаты. В любом случае, я уже думал о реализации адаптивного размера шага. Я попытаюсь решить это так. Еще раз спасибо за ответ.
Снимаю свой предыдущий комментарий. Я думаю, что мое внимание привлек автор «Последнее изменение», и я подумал, что Альфонсо написал вопрос и ответил на него. Мои извинения.
@ChrisR Не беспокойтесь!
+1Глядя на первый график ОП и график различий, я первым делом заметил, что синий всегда «нечеткий» по сравнению с красным; отклонения ошибок внутри орбиты кажутся большими. Это, безусловно, согласуется с тем, что размер шага слишком велик для используемого метода; усредненные по любой полной орбите ошибки могут иметь тенденцию к устранению, но внутриорбитальные симметричные отклонения могут сделать его более шумным.
@uhoh спасибо за ваш ответ. Возможно, «нечеткость» синего связана со спецификацией строки «--», установленной на графике Matlab. Несмотря на это, как бы вы решили? Я попробовал интегратор типа Рунге-Кутта-Фельберга, как рекомендовалось, но ничего не изменилось.
@Frank, мы можем заменить «нечеткий» на «толстый» или «широкий». Что ж, поскольку есть разногласия, мы можем знать, что одно или оба в значительной степени просто неверны. Я рекомендую изолировать интеграцию от всего остального. Получите начальный вектор состояния и запустите простую программу в Matlab или, что еще лучше, в Python с открытым исходным кодом, и интегрируйте одну, а затем несколько орбит и сравните напрямую. Я добавлю ответ о том, как это сделать через несколько часов.
В дополнение к комментарию @uhoh вы также можете попробовать проверить разницу между этим распространителем GMAT и ode45 Matlab (а также вашим распространителем). Я также обычно поощрял бы Python, но я думаю, что компания OP ищет это в Matlab.

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

Самая фундаментальная проверка работоспособности с помощью числового решателя: как изменится результат при изменении размера шага? Держите его как можно проще. Смоделируйте за один день и попробуйте разные размеры шагов. Мол, попробуйте шаг 6 секунд, 30 секунд, 200 секунд, 600 секунд. Чем они отличаются?

600 секунд наверное будут вести себя очень странно, а вот остальные - не знаю, надо попробовать. Если решение для h = 6 секунд по существу такое же, как ваше решение для h = 60 с, в частности, если оно постоянно имеет те же самые отклонения от GMAT, то у вас, очевидно, есть более фундаментальная проблема, не имеющая отношения к вообще интегратор. Это может быть много разного.

Только в том случае, если вы установили, что вы определенно настроили все так, как должно быть, вы должны сосредоточиться на том, чтобы получить максимальную отдачу от вашего численного решателя.

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

Ваш числовой интегратор РК4 не подходит для этой задачи. Порядок слишком низкий для требуемой точности. Это также фиксированный размер шага. Вы не провели исследование по этому вопросу. Скорость распространения ошибок для интеграторов RK с фиксированным шагом почти настолько плоха, насколько это возможно. Исследуйте интеграторы предикторов-корректоров высокого порядка... порядка 8 с (как минимум) 64-битными операциями с плавающей запятой (даже если они имеют фиксированный шаг).

Вы можете взглянуть на встроенные интеграторы RK, которые используют интегратор RK более низкого порядка в качестве предиктора в интеграторе RK более высокого порядка. Посмотрите на RK78 как на замену вашему RK4. В прошлом этот интегратор был рабочей лошадкой для орбитальных механиков. Это даст вам гораздо лучшее согласование данных и хорошую точность в течение более длительных интервалов времени.

Адамс-Моултон, или еще лучше, Адамс-Башфорт-Моултон, являются предикторами-корректорами, которые можно легко закодировать в любом порядке. Их скорость распространения ошибок меньше, чем у любого метода РК того же порядка. Интеграторы Гаусса-Радау также очень точны и хорошо подходят для интеграции орбитальной механики. Не кодируйте их сами. Скачать.

Существует множество методов предиктора-корректора с переменным шагом, которые обеспечат наилучшую точность и наименьшую скорость распространения ошибок в течение более длительных интервалов времени. С помощью этих методов ваши числовые интегрирования обеспечат наилучшую точность, которую вы можете ожидать... но они должны быть высокого порядка. Заказ 8 хорош.

Прежде чем приступить к проекту с любым интегратором, упомянутым здесь или где-либо еще, проверьте выбранный вами интегратор на решении задачи двух тел, выраженной в декартовых координатах. Это наихудший способ записи уравнений двух тел, и такая форма уравнений лучше всего нагрузит ваших интеграторов. Поскольку решение задачи двух тел имеет замкнутую форму в любой системе координат, вы увидите, как выбранные интеграторы справляются с распространением ошибок и точностью на коротких и длинных интервалах времени.

Ваш проект ни в коем случае не является тривиальным, и вам понадобится интегратор, который справится с этой задачей. Вы получите хорошее представление об интеграторе, который можно использовать для численного интегрирования задачи двух тел. Убедитесь, что ваши начальные условия дадут вам орбиту из двух тел с достаточно высоким эксцентриситетом: e >= 0,3. Не сходите с ума с очень высоким e (0,8). Интеграторы с фиксированным шагом быстро исчезнут с орбитами с большим эксцентриситетом (e = 0,8). Скорость их гибели обратно пропорциональна порядку интегратора.

Удачи.

Отличный ответ! Приятно видеть еще одного поклонника методов Адамса. Есть ли у вас какие-либо советы о том, когда выбрать Адамса-Башфорта-Моултона или Гаусса-Радау? О последних я слышал, но никогда ими не пользовался, а АБМ пользуюсь почти постоянно.
Спасибо за ваш ответ, я уже пробовал RK89 с фиксированным шагом, но результаты такие же. Также ничего не меняется с RK45 с адаптивным шагом. Но вы правы, мне нужно больше изучить методы этого интегратора и выбрать наиболее подходящий для решения проблемы. Чего я не могу объяснить, так это того, почему один и тот же распространитель хорошо работает в его форме Simulink, в то время как он показывает проблемы как скрипт Matlab.
Этот ответ, безусловно, содержит некоторые важные моменты, но они сформулированы слишком общим образом. Увеличение порядка решателя из-за того, что происходит что-то странное, чего я не понимаю, на самом деле довольно ужасная идея. Да, более высокий порядок, конечно, может дать лучшую точность, но он не решает фундаментальных проблем со стабильностью и уж точно не устраняет ошибки настройки. RK4 достаточно хорош, чтобы исследовать эти вещи, а более низкий порядок фактически делает его более подходящим для понимания того, как числовые интеграторы ведут себя с различными ODE и различными постоянными или адаптивными размерами шага.
Я согласен с @leftaroundabout на 1000%. Я бы сказал, что для таких несложных задач, как почти круговая орбита с квадрупольным членом всего 1 часть на тысячу, вы можете получить совершенно хорошие результаты с RK4 или даже ниже, если размер вашего шага подходит . (спасибо Дэвиду Хаммену ). Для такого рода задач интегратор более высокого порядка дает вам ту же точность с большим размером шага и, следовательно, потенциально более коротким временем вычислений, которое для этого будет составлять доли секунды или меньше.
@угу да. Хотя я бы не сказал, что «почти круговой» и «всего 1 ‰ квадруполь» являются причинами, по которым решатель низкого порядка подходит. Это скорее причины, по которым адаптивный размер шага ничего не даст. И это на самом деле также сценарий, в котором многошаговые методы высокого порядка могут действительно проявить себя, позволяя большой размер шага и низкие вычислительные затраты для вычисления миллионов орбит в будущем. ... только это довольно бессмысленно, если результаты уже являются фиктивными после всего лишь нескольких орбит ... Для вычислений в будущем в менее регулярных ситуациях может быть лучше симплектический метод более низкого порядка.
@leftaroundabout оба могут быть правдой, они не являются исключительными; но да, конечно, адаптивный размер шага в этом случае мало помогает. Если бы эксцентриситет был 0,9 или величина квадрупольного ускорения была бы аналогична монопольному, а размер шага был фиксированным 30 секунд, RK4 был бы паршивым, и экзотический более высокий порядок справился бы лучше, но да, вы правы, RK45 подкопался бы и понизил бы это размер шага там, где это необходимо. Сейчас работаю над ответным постом...
Здесь мы снова имеем мнения, а не факты о численном интегрировании траекторий орбитальной механики. Я полагаю, что ни у кого нет опыта распространения ошибок различных методов интегрирования. RK4 НИКОГДА не используется в каком-либо серьезном численном исследовании проблемы орбитального меха. Распространение ошибок усечения является основным виновником. То, о чем вы хотите беспокоиться, - это распространение ошибки округления, и если это вас беспокоит, это лучшее (наименьшее беспокойство), которое вы можете получить.
Кроме того, e = 0,9 не представляет реального или академического интереса. Размер шага в 30 секунд с радостью уничтожит ваш вывод RK4 благодаря быстрому накоплению ошибок усечения. Основная проблема с интеграторами низкого порядка — распространение ошибки усечения. Ошибка округления никогда не будет проблемой. Вы хотите использовать интегратор, если вас беспокоит только распространение ошибки округления. Зачем использовать 64-битные числа с плавающей запятой и интегратор младшего порядка? Небольшие размеры шага для компенсации интеграторов низкого порядка с фиксированным шагом никогда не являются приемлемым решением из-за быстрого нарастания ошибки усечения.
Некоторые отличные интеграторы с переменным шагом для задач об орбитах используют скорость изменения кривизны трехмерной траектории для определения изменений размера шага. Добавление таких соображений в код интегратора включает в себя очень несколько дополнительных строк кода с использованием значений переменных состояния, полученных на предыдущем шаге. Переменный шаг помогает свести к минимуму распространение ошибок. Почему бы не использовать его? Это не требует больших затрат времени на выполнение.
Райан К.: Мне лично не нравится метод Гаусса-Радау из-за его сложности. Это сложный код (по крайней мере, для меня). Мне нравится добавлять функции в готовый код, а это сложно сделать с помощью Гаусса-Радау. ABM легко написать самостоятельно. Вы можете создать код ABM с входными данными, определяющими порядок метода. Я экспериментировал с этим сквозным порядком 25 (слишком много!) и 128-битной плавающей запятой. ABM очень эффективен, и вы можете настроить порядок так, чтобы с каждым шагом учитывать ошибку округления машины. С дополнительной настройкой вы можете изменить порядок во время интеграции по ходу дела.
Плюс один только за «Не кодируйте их сами. Загрузите».
@AngePurs Проблема с ABM заключается в том, что он интегрирует ODE первого порядка. Хотя ОДУ второго порядка можно преобразовать в первый порядок, при этом геометрия отбрасывается. Существуют варианты ABM, которые рассматривают положение и скорость отдельно, такие как Störmer-Cowell и Gauss-Jackson. Примечание: слишком высокий порядок в интеграторе (не путать с порядком ОДУ) может привести к потере точности и стабильности. 25-й порядок слишком высок, если только вы не используете арифметику с чрезвычайно высокой точностью - и тогда вы платите высокую вычислительную цену за использование арифметики с чрезвычайно повышенной точностью.
Различия между числовыми пропагаторами
СпросиСетьПолитика конфиденциальности1y, 0v