12 Mar 2025
Есть целая куча методов численного интегрирования в самыми разными свойствами. Ниже примеры и краткий обзор алгоритмов.
Свойства численных методов
-
Порядок точности. Описывает то, как оценивать ошибку интегрирования. Первый порядок точности - ошибка зависит от второй производной, второй порядок точности - от третьей. Методы первого порядка точности лучше вообще не использовать, большинство более-менее приличных методов имеют второй порядок точности, самые замороченные - четвёртый. Теоретически можно больше, но я не пробовал.
-
Симметричность относительно времени. Например, в интегрировании Верле можно сделать несколько шагов “вперёд”, потом столько же “назад” и оказаться ровно в той же точке.
-
Симплектические интеграторы - звучит страшно, но суть в том, что при моделировании тел будет сохраняться энергия системы.
-
Одно и много-шаговые методы. В многошаговых один “шаг” физики требует вычисления сил в нескольких точках.
-
Явные, неявные и полу-неявные методы.
- Явные методы - самые простые, где следующая точка считается на основе предыдущих.
- Неявные - ответ в уравнении есть и слева и справа, например y(t+dt) = y(t) + y’(t+dt) dt. В общем случае решаются итеративно - раз за разом находим y(t+dt) и подставляем это в правую часть уравнения. Обычно более устойчивые, чем явные методы.
- Полу-явные - интересный компромисс, сочетающий в себе простоту явного метода и при этом лучшую сходимость.
Численные методы
Метод Стрёмера-Верле
Простой и очень интересый метод. Второй порядок точности, одношаговый, симплектический и симметричный.
Скорость в явном виде не хранится, вместо неё - хранятся текущая и предыдущая позиция.
\[y(t+dt) = 2 y(t) - y(t - dt) + y''(t) \dfrac{dt^2}{2}\]Шаг времени динамически поменять не получится, хотя мне кажется, что при некоторых ухищрениях можно.
Например, найти скорость в точке:
\[y'(t) = \dfrac{y(t+dt) - y(t-dt)}{2dt}\]И потом каким-нибудь другим точным методом интегрирования на основе y(t) и y’(t) получить следующую точку через произвольный интервал.
Ещё интересный момент - если есть ограничения типа “ниточек” или “палочек” между точками, то можно просто брать и двигать обе точки обратно пропорционально их массам, тогда и ограничение удовлетворится и точность не пострадает. Скорости в явном виде нет, а пересчёт координат будет очень простой.
Полу-явный метод Эйлера
Неточный, но зато является хорошей иллюстрацией полу-явного метода. Симплектический и сходится лучше, чем явный метод Эйлера. Явный метод Эйлера - не симплектический, и его вообще нет смысла использовать.
Допустим, что у нас в системе есть координаты и скорости. Вместо того, чтобы обновить сразу и координату и скорость, обновляем сначала что-то одно, а при вычислении второго используем более новую версию.
\[x_{n+1} = x_n + f(t_n, v_n) dt\] \[v_{n+1} = v_n + g(t_n, x_{n+1}) dt\]Можно сделать любым способом - либо считать новую скорость сразу для новой позиции, либо новую позицию на основе новой скорости. Вычислений столько же, сколько в явном методе Эйлера, но из-за симплектичности точность при симуляции тел выше. Но порядок всё ещё первый.
Метод средней точки
Явный метод второго порядка.
Идея простая - делаем “подшаг” на половину времени вперёд, считаем силы там, и уже на основе этой промежуточной силы двигает тело из начальной точки в следующую.
\[y_{n+1} = y_{n} + f(t_n + \dfrac{dt}{2}, y_n + \dfrac{dt}{2}f(t_n, y_n)) dt\]Есть неявная вариация метода:
\[y_{n+1} = y_n + f(t + \dfrac{dt}{2}, \dfrac{y_n + y_{n+1}}{2}) dt\]Здесь надо сделать несколько итераций. Средняя точка выражается как “середина” между начальной и конечной, а ускорение в средней точке должно задать перемещение из начальной в конечную.
Heun’s method
Похож на метод средней точки, тоже явный и второго порядка.
Сдесь используется чуть другая идея - делается шаг вперёд на dt, там считается новая сила, а дальше снова делается шаг от начальной точки, но уже с усреднённой силой в начале и в конце.
\[y_{n+1} = y_n + \dfrac{f(y_n) + f(y_n + f(y_n)dt)}{2} dt\]В принципе, можно повторять шаги несколько раз и метод станет немножко похож на неявный метод средней точки.
Классический явный метод Рунге-Кутты четвёртого порядка.
Вообще очень многие методы можно назвать методами Рунге-Кутты, это такое обобщение для всего. Вот список
Требует вычисленя сразу в чётырёх точках, даёт четвёртый порядок точности.
Классический метод на мой взгляд самый удобный
\[y_{d1}' = f(y_n, t_n)\] \[y_{d2}' = f(y_n + \dfrac{dt}{2} y_{d1}', t_n + \dfrac{dt}{2})\] \[y_{d3}' = f(y_n + \dfrac{dt}{2} y_{d2}', t_n + \dfrac{dt}{2})\] \[y_{d4}' = f(y_n + \dfrac{dt} y_{d3}', t_n + dt)\] \[y_{n+1} = y_n + \dfrac{y_{d1}' + 2 y_{d2}' + 2 y_{d3}' + y_{dt}'}{6} dt\]Есть вариация, где используются другие коэффициенты типа 1/3, -1, 1, 3/8, но я не увидел в нём большого смысла, так как вычисления выглядят чуть сложнее, а точность почти такая же.
Неявный метод Гаусса-Лежандра четвёртого порядка
Отличается от предыдущего метода тем, что он более устойчивый. Количество итераций любое. У меня получилось, что на трёх итерациях точность как у Рунге-Кутты, но вот если явный метод не справляется из-за слишком большого шага времени, то этим методом можно сделать побольше итераций и вытянуть точность.
Я не возьмусь писать формулы, но можно почитать википедию или посмотреть мой код.
На мой взгляд чуда не случается. В каких-то случаях может быть лучше сделать штуки 4 итерации, чем уменьшать шаг рассчёта в два раза и считать явным методом, но разница не очень большая. Наибольший эффект даёт переход от метода второго порядка точности к любому методу четвёртого порядка
Итоги
Методы первого порядка точности использовать не надо, точность ужасна.
Если второго порядка точности хватает, я бы советовал метод Верле из-за обилия его интересных особенностей. Либо какой-либо из перечисленных выше методов второго порядка. Они все примерно одинаковые по сложности и точности.
Если нужно точнее, то лучше сразу брать явный метод Рунге-Кутты четвёртого порядка либо неявный метод Гаусса-Лежандра.
Я проверял методы на том, как они симулируют вращение тел с прецессией. код тут. Рядом в папке test можно найти тесты к ним.