Мне было просто любопытно, почему методы Рунге-Кутты высокого порядка (т.е. больше 4) почти никогда не обсуждаются / не используются (по крайней мере, насколько мне известно). Я понимаю, что это требует большего вычислительного времени на шаг (например, RK14 со встроенным шагом 12-го порядка ), но есть ли другие недостатки использования методов Рунге-Кутты более высокого порядка (например, проблемы стабильности)? Применительно к уравнениям с сильно колеблющимися решениями в экстремальных временных масштабах, не будут ли предпочтительны такие методы более высокого порядка?

Существуют тысячи статей и сотни кодов, использующих методы Рунге-Кутты пятого порядка или выше. Обратите внимание, что наиболее часто используемый явный интегратор в MATLAB - это ODE45, который продвигает решение, используя метод Рунге-Кутты 5-го порядка.

Примеры широко используемых методов Рунге-Кутты высокого порядка

В статье Dormand & Prince, в которой описан метод 5-го порядка, содержится более 1700 ссылок, согласно Google Scholar . Большинство из них - документы, использующие их метод для решения какой-то проблемы. У метода Cash-Karp более 400 ссылок . Возможно, наиболее широко используемый метод порядка выше 5 - это метод Принца-Дормана 8-го порядка, который содержит более 400 ссылок в Google Scholar . Я мог бы привести много других примеров; и имейте в виду, что многие (если не большинство) людей, использующих эти методы, никогда не цитируют статьи.

Отметим также, что методы экстраполяции и отложенной коррекции высокого порядка являются методами Рунге-Кутты .

Методы высокого порядка и ошибка округления

Если ваша точность ограничена ошибками округления, вам следует использовать метод более высокого порядка . Это связано с тем, что методы более высокого порядка требуют меньше шагов (и меньше оценок функций, даже если на шаг приходится больше оценок), поэтому они допускают меньше ошибок округления. Вы можете легко проверить это самостоятельно с помощью простых экспериментов; это хорошая домашняя задача для первого курса численного анализа.

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

Методы высокого порядка могут быть такими же стабильными

$A$$B$

Методы высокого порядка в небесной механике

Ты спрашиваешь

Применительно к уравнениям с сильно колеблющимися решениями в экстремальных временных масштабах, не будут ли предпочтительны такие методы более высокого порядка?

Вы совершенно правы! Ярким примером этого является небесная механика. Я не эксперт в этой области. Но эта статья , например, сравнивает методы небесной механики и даже не рассматривает порядок ниже 5. В ней делается вывод, что методы порядка 11 или 12 часто являются наиболее эффективными (с методом Принца-Дормана порядка 8 также часто очень эффективная).

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

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

Просто чтобы добавить к отличному ответу Брайана Борчера, многие реальные приложения допускают очень жесткие ODE или DAE. Интуитивно понятно, что эти проблемы со временем претерпевают негладкие, резкие изменения, поэтому их лучше моделировать с помощью полиномов низкого порядка, которые распределены по небольшим размерам шагов, в отличие от полиномов высокого порядка, растянутых по длинным шагам. Кроме того, стабильность часто требует использования неявных методов, для которых вычислительные потери методов более высокого порядка намного круче.

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

$r\le2$ . Среди всех многошаговых методов порядка 2 правило трапеции имеет наименьшую постоянную ошибки.

Аналогичные (но гораздо более сложные) утверждения можно сделать для L-устойчивости в формулах РК. Во всех случаях увеличение порядка часто не всегда приводит к более точным решениям. Ниже приводится выдержка из оригинальной статьи Протеро и Робинсона 1974 года:

Используя A-стабильные одношаговые методы для решения больших систем жестких нелинейных дифференциальных уравнений, мы обнаружили, что
(a) некоторые A-стабильные методы дают весьма нестабильные решения и
(b) точность решений, полученных, когда уравнения Жесткость часто оказывается не связанной с порядком используемого метода.

Более подробное описание этой темы см. В классическом тексте Хайрера и Ваннера «Решение обыкновенных дифференциальных уравнений II: жесткие и дифференциально-алгебраические задачи», 1991 г.

На практике жесткие уравнения почти всегда решаются с использованием правила трапеции или формулы TR-BDF2 (функции ode23t и ode23tb в MATLAB). Оба из них являются неявными методами второго порядка. Конечно, там, где стабильность не является проблемой (то есть в нестационарных уравнениях), мы можем выбирать из нескольких вариантов; RK45 - самый распространенный выбор.

Настройка Benchmark

В программном обеспечении Julia DifferentialEquations.jl мы реализовали множество методов более высокого порядка, включая методы Феагина. Вы можете увидеть это в нашем списке методов , а затем есть множество других, которые вы можете использовать в качестве поставляемых таблиц . Поскольку все эти методы собраны вместе, вы можете легко сравнить их между собой. Вы можете увидеть тесты, которые у меня есть в сети , и увидеть, что очень просто тестировать множество различных алгоритмов. Так что если вы хотите потратить несколько минут на тесты, сделайте это. Вот краткое изложение того, что выходит.

Прежде всего, важно отметить, что, если вы посмотрите на каждый из тестов, вы увидите, что наши DP5(Порядок Дормана-Принса 5) и DP8методы работают быстрее, чем коды ( dopri5и dop853) Hairer Fortran , и поэтому эти реализации очень хорошо оптимизированы. , Они показывают, что, как отмечалось в другом потоке, чрезмерное использование методов Dormand-Prince происходит потому, что методы уже написаны, а не потому, что они все еще являются лучшими. Таким образом, реальное сравнение между наиболее оптимизированными реализациями проводится между методами Tsitorous, методами Verner и Feagin из DifferentialEquations.jl.

Результаты

В целом, методы порядка выше 7 имеют дополнительные вычислительные затраты, которые обычно не перевешиваются по порядку, учитывая выбранные допуски. Одна из причин этого заключается в том, что выбор коэффициентов для методов более низкого порядка более оптимизирован (они имеют небольшие «основные коэффициенты ошибки усечения», которые имеют большее значение, когда вы не асимтотически малы). Вы можете видеть, что во многих задачах, подобных здесь, методы Verner Efficient 6 и 7 работают очень хорошо, но такие методы, как Verner Efficient 8, могут иметь меньший наклон. Это связано с тем, что «выгоды» более высокого порядка составляют более низкие допуски, поэтому всегда существует допуск, где методы более высокого порядка будут более эффективными.

Однако тогда возникает вопрос, насколько низко? В хорошо оптимизированной реализации это становится довольно низким по двум причинам. Первая причина в том, что методы низшего порядка реализуют нечто, называемое FSAL (первый такой же, как последний). Это свойство означает, что методы более низкого порядка повторно используют оценку функции из предыдущего шага на следующем шаге и, таким образом, фактически имеют на одну меньше оценку функции. Если это используется правильно, то что-то вроде метода 5-го порядка (Tsitorous или Dormand-Prince) фактически берет 5 оценок функций вместо 6, которые предложили бы таблицы. Это также верно для метода Вернера 6.

Другая причина связана с интерполяцией. Одна из причин использования метода очень высокого порядка состоит в том, чтобы делать меньше шагов и просто интерполировать промежуточные значения. Однако для получения промежуточных значений интерполяционной функции может потребоваться больше оценок функции, чем было использовано для выполнения шага.Если вы посмотрите на методы Вернератребуется 8 дополнительных функций для метода Порядка 8, чтобы получить интерполант Порядка 8. Во многих случаях методы низкого порядка предоставляют «свободный» интерполант, например, большинство методов 5-го порядка имеют свободную интерполяцию 4-го порядка (без дополнительных вычислений функции). Таким образом, это означает, что если вам нужны промежуточные значения (которые вам понадобятся для хорошего графика, если вы используете метод высокого порядка), есть некоторые дополнительные скрытые затраты. Фактор в том, что эти интерполированные значения действительно важны для обработки событий и решения дифференциальных уравнений с задержкой, и вы понимаете, почему дополнительные затраты на интерполяцию влияют на.

Так что насчет методов Феагина?

Таким образом, вы увидите, что методы Feagin подозрительно отсутствуют в тестах. Они хороши, тесты сходимости работают с произвольными точными числами и т. Д., Но для того, чтобы заставить их работать хорошо, нужно попросить довольно абсурдно низкие допуски. Например, я обнаружил в неопубликованных тестах, что они Feagin14превосходят Vern9(эффективный метод Вернера 9-го порядка) с такими допусками, как 1e-30. Для приложений с хаотической динамикой (например, в задачах Плайдеса или астрофизики с тремя телами) вам может потребоваться такой уровень точности из-за чувствительной зависимости (ошибки в хаотических системах быстро увеличиваются). Тем не менее, большинство людей, вероятно, работают с числами с плавающей запятой двойной точности, и я не нашел эталон, где они превзошли бы в этой области терпимости.

Кроме того, нет никакого интерполента, чтобы согласиться с методами Феагина. Так что я просто делаю на них интерполяцию Эрмита третьего порядка так, чтобы она существовала (и она работает на удивление хорошо). Однако, если нет стандартной интерполяционной функции, вы можете выполнить рекурсивный метод Эрмита (используйте эту интерполяцию для получения средней точки, затем выполните интерполяцию 5-го порядка и т. Д.), Чтобы получить интерполяцию высокого порядка, но это очень дорого, и в результате Интерполяция не обязательно имеет низкий коэффициент ошибки усечения принципа (так что это хорошо, только когда dtдействительно мало, что является полной противоположностью того случая, который мы хотим!). Поэтому, если вам когда-нибудь понадобится действительно хорошая интерполяция, чтобы соответствовать вашей точности, вам нужно, по крайней мере, вернуться к чему-то подобному Vern9.

Примечание об экстраполяции

Обратите внимание, что методы экстраполяции - это просто алгоритмы для генерации методов Рунге-Кутты произвольного порядка. Однако для их порядка они делают больше шагов, чем необходимо, и имеют высокие коэффициенты ошибки усечения, и поэтому они не так эффективны, как хорошо оптимизированный метод RK в данном порядке. Но, учитывая предыдущий анализ, это означает, что существует область крайне низкой толерантности, где эти методы будут работать лучше, чем «известные» методы РК. Но в каждом тесте, который я проводил, кажется, что я не достиг такого низкого уровня.

Примечание о стабильности

Выбор действительно не имеет ничего общего с проблемами стабильности. На самом деле, если вы пройдете через таблицу DifferentialEquations.jl (вы можете только plot(tab)для областей стабильности), вы увидите, что большинство методов имеют подозрительно похожие области стабильности. Это на самом деле выбор. Обычно при выводе методов автор обычно делает следующее:

1. Найти самые низкие основные коэффициенты ошибки усечения (то есть коэффициенты для членов следующего порядка)
2. С учетом ограничений порядка
3. И приблизить область стабильности к методу Дормана-Принца 5-го порядка.

Почему последнее условие? Что ж, потому что этот метод всегда стабилен в зависимости от того, как выполняется адаптивный выбор размера шага, управляемый PI, так что это хорошая планка для «достаточно хороших» областей стабильности. Так что не случайно все регионы стабильности имеют тенденцию быть похожими.

Вывод

Есть компромиссы в каждом выборе метода. Методы RK самого высокого порядка просто не так эффективны при более низких допусках, потому что оптимизировать выбор коэффициентов труднее, а также потому, что число оценок функций составляет (и растет даже быстрее, когда задействованы интерполяции). Однако, если допуск становится достаточно низким, они выигрывают, но требуемые допуски могут быть намного ниже «стандартных» применений (то есть действительно применимы только к хаотическим системам).