Ищем Рунге-Кутта 8-го порядка в C / C ++

10

Я хотел бы использовать метод 8-го порядка Рунге-Кутты (89) в приложении по небесной механике / астродинамике, написанном на C ++, на машине Windows. Поэтому мне интересно, знает ли кто-нибудь хорошую библиотеку / реализацию, которая документирована и бесплатна для использования? Это нормально, если он написан на C, если нет никаких проблем с компиляцией.

До сих пор я нашел эту библиотеку (mymathlib) . Код кажется нормальным, но я не нашел никакой информации о лицензировании.

Можете ли вы помочь мне, раскрывая некоторые альтернативы, которые вы, возможно, знаете и которые подойдут для моей проблемы?

РЕДАКТИРОВАТЬ:
я вижу, что на самом деле не так много исходных кодов C / C ++, как я ожидал. Поэтому версия Matlab / Octave тоже подойдет (все еще должна быть бесплатной для использования).

Джеймс С
источник

Ответы:

8

И GNU Scientific Library (GSL) (C) и увеличить Odeint (C ++) функция 8 - го порядка методов Рунге-Кутта.

Оба с открытым исходным кодом, а в Linux и Mac они должны быть непосредственно доступны из менеджера пакетов. Под окнами вам, вероятно, будет проще использовать Boost, а не GSL.

GSL публикуется под лицензией GPL, а Boost Odeint под лицензией Boost.

Изменить: Хорошо, Boost Odeint не имеет метод Runge-Kutta 89, только 78, но он обеспечивает рецепт для создания произвольных степперов Runge-Kutta.

Однако методы 8-го порядка довольно высоки и, скорее всего, излишни для вашей проблемы.

Prince-Dormand относится к определенному виду Runge-Kutta и не имеет прямого отношения к порядку, но наиболее распространенным является 45. Matlabs ode45, рекомендуемый алгоритм ODE, является реализацией Prince-Dormand 45. Это тот же алгоритм, который реализован в Boost Odeint Runge_Kutta_Dopri5 .

LKlevin
источник
1
Спасибо за ответ. Хорошо, что сейчас стыдно, я взглянул на Boost Odeint еще до того, как спрашивать здесь, и нашел только «runge_kutta_fehlberg78». Это правильно? На самом деле я не знаю различий между метондами, когда они используются на практике, но я искал RK89 (он также называется Dormand-Prince, когда я ищу в Интернете). Можете ли вы прокомментировать или расширить свой ответ по этому вопросу, пожалуйста? Спасибо.
Джеймс C
Обновленный пост, чтобы ответить на ваши вопросы. Prince-Dormand 45, скорее всего, хорошо решит ваши проблемы.
Л.Клевин
15

Если вы занимаетесь небесной механикой в ​​течение длительного времени, использование классического интегратора Рунге-Кутты не сохранит энергию. В этом случае лучше использовать симплектический интегратор. Boost.odeint также реализует симплектическую схему Рунге-Кутты 4-го порядка, которая будет работать лучше в течение длительных интервалов времени. Насколько я могу судить, GSL не реализует никаких симплектических методов.

Джефф Оксберри
источник
Спасибо за ответ. Может ли симплектическая Рунге-Кутта 4-го порядка дать лучшие результаты, чем RKF78, если она используется со спутниками Земли (как на низкой, так и на более глубокой космической орбите), возможно, в течение периода 1-3 орбит?
Джеймс С.
@JamesC Да. В течение длительного периода симплектический метод намного лучше.
eccstartup
@eccstartup - Что бы вы посчитали здесь долгим периодом? Потому что это может быть одна орбита планеты вокруг Солнца или несколько орбит метеорологического спутника вокруг Земли и т. Д.
Джеймс С.
@JamesC Я не заметил этой большой проблемы. Но для моих модельных задач с вычислением многих орбит симплектические методы дают очень совершенные орбиты.
eccstartup
Итак, это совет для программирования вашей собственной версии неявного метода Рунге-Кутты, который включает в себя множество симплектических методов с таким высоким порядком, как вам нужно.
eccstartup
4

резюмируя некоторые моменты:

  1. Если это долгосрочная интеграция недиссапативной модели, то вам нужен симплектический интегратор.
  2. В противном случае, поскольку это уравнение движения, методы Рунге-Кутты Нистрома будут более эффективными, чем преобразование в систему первого порядка. Существуют методы RKN высокого порядка благодаря DP. Есть несколько реализаций, как здесь, в Джулии, они документированы, а вот MATLAB .
  3. Методы Рунге-Кутты высокого порядка необходимы только в том случае, если вам нужно высокоточное решение. Если это более низкие допуски, то РК 5-го порядка, вероятно, будет быстрее (для той же ошибки). Лучшее, что нужно сделать, если вам часто приходится решать эту проблему, - это протестировать несколько различных методов. В этом наборе тестов для задач из трех тел мы видим, что (для той же ошибки) методы RK высокого порядка являются лишь незначительным улучшением скорости, хотя, как ошибка -> 0, вы можете видеть, что улучшение уже достигает> 5x по сравнению с Dormand - Принять 45 ( DP5), когда вы смотрите на 4 цифры точности (хотя для этого допуски значительно ниже. Допуски являются лишь приблизительным критерием в любой проблеме). По мере того, как вы увеличиваете допуски, даже улучшение улучшается от метода РК высокого порядка, но вам, возможно, придется начать использовать числа с более высокой точностью.
  4. Алгоритм порядка Дормана-Принса 7/8 имеет другую схему 8-го порядка, чем метод DP853 в Hairer dop853и DifferentialEquations.jl DP8(которые совпадают). Последний метод 853 не может быть реализован в стандартной табличной версии метода Рунге-Кутты, поскольку его оценка ошибок является нестандартной. Но этот метод гораздо более эффективен, и я бы не рекомендовал даже использовать более старые методы Fehlberg 7/8 или DP 7/8.
  5. Для методов высокого порядка РК методы Вернера "Эффективные" являются золотым стандартом. Это проявляется в тестах, которые я связал. Вы можете закодировать их в Boost самостоятельно или использовать один из двух пакетов, которые их реализуют, если вам это нужно (Mathematica или DifferentialEquations.jl).
Крис Ракауцкас
источник
2

Я хотел бы добавить, что хотя то, что Джефф Оксберри предлагает для долгосрочной интеграции (с использованием симплектических интеграторов), верно, в некоторых случаях это не сработает. Более конкретно, если у вас есть диссипативные силы, ваша система больше не сохраняет энергию, и поэтому вы не можете прибегнуть к симплектическим интеграторам в этом случае. Человек, задающий вопрос, говорил о низких околоземных орбитах, и такие орбиты демонстрируют большое атмосферное сопротивление, то есть диссипативную силу, которая препятствует использованию таких симплектических интеграторов.

В этом конкретном случае (а также в тех случаях, когда вы не можете использовать / не имеете доступа / не хотите использовать симплектические интеграторы), я бы рекомендовал использовать интегратор Bulirsch-Stoer, если вам нужна точность и эффективность в течение длительного периода времени. Он хорошо работает по опыту, а также рекомендуется в Числовых рецептах (Press et al., 2007).

Viiv
источник
Нет, не рекомендую Численные Рецепты. Особенно в большинстве случаев не рекомендуется рекомендовать Burlirsch-Stoer. Это известная проблема с книгой. Смотрите кучу опровержений от ведущих исследователей в этой области здесь: uwyo.edu/buerkle/misc/wnotnr.html . Если вам нужны тесты по этой теме, посмотрите первую книгу Hairer, где вы увидите, что BS почти никогда не преуспевает. Более высокий порядок эффективнее только тогда, когда ошибки достаточно малы, и мы (и другие) провели сравнительный анализ, чтобы достаточно последовательно показать, что он эффективен только для точности с плавающей запятой.
Крис
Я не могу говорить за NR слишком много, так как я использовал его в основном для ODE, но мне кажется, что жалобы на странице, на которую вы ссылаетесь, старые и были отмечены авторами NR в своем ответе (конец страницы), но это не по теме. Что касается долгосрочной интеграции орбит с высокой точностью (скажем, 13-14 цифр), о которой я упоминал в своем ответе, то уже давно доказано, что методы экстраполяции работают хорошо (см. Главу Монтенбрюка и Джилла о числовой интеграции). Более поздние статьи также используют его, и он доказал мне и другим надежный и эффективный метод.
Viiv
M & G тестирует его только на dop853, и более современные методы высокого порядка RK, такие как благодаря Вернеру, гораздо эффективнее. M & G также, кажется, измеряет только с помощью оценки функций, которые являются слабым индикатором времени. Это также не время против методов Рунге-Кутты Нистрома, которые специально предназначены для ODE 2-го порядка и более эффективны, чем методы RK первого порядка, применяемые ко второму порядку. В 13-14 разрядах BS, вероятно, конкурентоспособна по большинству задач, но это далеко от очевидного выбора, и я не видел диаграммы точности работы с недавними методами, не согласными с этим.
Крис
M & G проводят тестирование RKN против RK, а BS и других - против RKN (стр. 123-132 и 151-154) и говорят, что они являются наиболее эффективными из методов RK (не считая Вернера, даже если они цитируют его). Показано, что BS эффективна при 13-14 цифрах, что было моим утверждением, я видел, что она проверена на dop853, ABM (12), Taylor и стандарте RK8, и она хорошо работает. Я должен признать, что я не видел, чтобы он снова тестировался в RKN, но, как я вижу из M & G, он находится недалеко от FILG11, например. Я искренне заинтересован в РК Вернера и посмотрю ваши ссылки выше. У вас есть бумага, которая проверяет их все, чтобы увидеть?
Viiv
Я вернулся и повторно выполнил несколько тестов в DiffEqBenchmarks.jl и odexне очень хорошо справляется . Так что, по крайней мере, для ODE 1-го порядка и для допусков >=1e-13экстраполяция выглядит не очень хорошо, и обычно она даже не близка. Это соответствует требованию выше.
Крис