Интеграция Верлет

Интеграция Верлета ( Французское произношение: [vɛʁˈlɛ] ) — численный метод, используемый для интегрирования Ньютона уравнений движения . Он часто используется для расчета траекторий частиц в молекулярно-динамическом моделировании и компьютерной графике . Алгоритм был впервые использован в 1791 году Жаном Батистом Деламбром и с тех пор много раз открывался заново, последний раз Лу Верле в 1960-х годах для использования в молекулярной динамике. Он также использовался П. Х. Коуэллом и А.К.К. Кроммелином в 1909 году для расчета орбиты кометы Галлея и Карлом Стормером в 1907 году для изучения траекторий электрических частиц в магнитном поле (поэтому его также называют методом Штёрмера ). Интегратор Верле обеспечивает хорошую численную стабильность , а также другие свойства, которые важны в физических системах, такие как обратимость во времени и сохранение симплектической формы в фазовом пространстве , без значительных дополнительных вычислительных затрат по сравнению с простым методом Эйлера .

Базовый Штёрмер-Верлет

Для дифференциального уравнения второго порядка типа ${\ddot {\mathbf {x} }}(t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}$ с начальными условиями $\mathbf {x} (t_{0})=\mathbf {x} _{0}$ и ${\dot {\mathbf {x} }}(t_{0})=\mathbf {v} _{0}$, приближенное численное решение $\mathbf {x} _{n}\approx \mathbf {x} (t_{n})$ в то время $t_{n}=t_{0}+n\,\Delta t$ с размером шага $\Delta t>0$ можно получить следующим методом:

1. набор ${\textstyle \mathbf {x} _{1}=\mathbf {x} _{0}+\mathbf {v} _{0}\,\Delta t+{\tfrac {1}{2}}\mathbf {A} (\mathbf {x} _{0})\,\Delta t^{2}}$,
2. для n = 1, 2,... итерация $$\mathbf {x} _{n+1}=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {A} (\mathbf {x} _{n})\,\Delta t^{2}.$$

Уравнения движения

Уравнение движения Ньютона для консервативных физических систем имеет вид

${\boldsymbol {M}}{\ddot {\mathbf {x} }}(t)=F{\bigl (}\mathbf {x} (t){\bigr )}=-\nabla V{\bigl (}\mathbf {x} (t){\bigr )},$

или индивидуально

$m_{k}{\ddot {\mathbf {x} }}_{k}(t)=F_{k}{\bigl (}\mathbf {x} (t){\bigr )}=-\nabla _{\mathbf {x} _{k}}V\left(\mathbf {x} (t)\right),$

где

• $t$ это время,
• $\mathbf {x} (t)={\bigl (}\mathbf {x} _{1}(t),\ldots ,\mathbf {x} _{N}(t){\bigr )}$ представляет собой ансамбль вектора положения $N$ объекты,
• $V$ скалярная потенциальная функция,
• $F$ – отрицательный градиент потенциала , дающий ансамбль сил, действующих на частицы,
• ${\boldsymbol {M}}$ — массовая матрица , обычно диагональная с блоками с массой $m_{k}$ для каждой частицы.

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

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

${\ddot {\mathbf {x} }}(t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}$

с некоторой подходящей векторной функцией $\mathbf {A} (\mathbf {x} )$ представляющее ускорение, зависящее от положения. Обычно исходное положение $\mathbf {x} (0)=\mathbf {x} _{0}$ и начальная скорость $\mathbf {v} (0)={\dot {\mathbf {x} }}(0)=\mathbf {v} _{0}$ также даны.

Интеграция Верлета (без скоростей)

Чтобы дискретизировать и численно решить эту задачу начального значения , шаг по времени $\Delta t>0$ выбирается, и последовательность точек отбора проб $t_{n}=n\,\Delta t$ обдуманный. Задача — построить последовательность точек $\mathbf {x} _{n}$ которые точно соответствуют пунктам $\mathbf {x} (t_{n})$ на траектории точного решения.

Там, где метод Эйлера использует аппроксимацию прямой разности первой производной в дифференциальных уравнениях первого порядка, интеграцию Верле можно рассматривать как использование аппроксимации центральной разности второй производной:

${\begin{aligned}{\frac {\Delta ^{2}\mathbf {x} _{n}}{\Delta t^{2}}}&={\frac {{\frac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}-{\frac {\mathbf {x} _{n}-\mathbf {x} _{n-1}}{\Delta t}}}{\Delta t}}\\[6pt]&={\frac {\mathbf {x} _{n+1}-2\mathbf {x} _{n}+\mathbf {x} _{n-1}}{\Delta t^{2}}}=\mathbf {a} _{n}=\mathbf {A} (\mathbf {x} _{n}).\end{aligned}}$

Интегрирование Верле в форме, используемой в качестве метода Штёрмера. использует это уравнение для получения следующего вектора положения из двух предыдущих без использования скорости как

${\begin{aligned}\mathbf {x} _{n+1}&=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {a} _{n}\,\Delta t^{2},\\[6pt]\mathbf {a} _{n}&=\mathbf {A} (\mathbf {x} _{n}).\end{aligned}}$

Ошибка дискретизации

Симметрия времени, присущая методу, снижает уровень локальных ошибок, вносимых в интегрирование за счет дискретизации за счет удаления всех членов нечетной степени, в данном случае членов в $\Delta t$ третьей степени. Локальная ошибка определяется количественно путем вставки точных значений. $\mathbf {x} (t_{n-1}),\mathbf {x} (t_{n}),\mathbf {x} (t_{n+1})$ в итерацию и вычисление разложений Тейлора за раз $t=t_{n}$ вектора положения $\mathbf {x} (t\pm \Delta t)$ в разных временных направлениях:

${\begin{aligned}\mathbf {x} (t+\Delta t)&=\mathbf {x} (t)+\mathbf {v} (t)\Delta t+{\frac {\mathbf {a} (t)\Delta t^{2}}{2}}+{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}\left(\Delta t^{4}\right)\\\mathbf {x} (t-\Delta t)&=\mathbf {x} (t)-\mathbf {v} (t)\Delta t+{\frac {\mathbf {a} (t)\Delta t^{2}}{2}}-{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}\left(\Delta t^{4}\right),\end{aligned}}$

где $\mathbf {x}$ это позиция, $\mathbf {v} ={\dot {\mathbf {x} }}$ скорость, $\mathbf {a} ={\ddot {\mathbf {x} }}$ ускорение и $\mathbf {b} ={\dot {\mathbf {a} }}={\overset {\dots }{\mathbf {x} }}$ рывок . (третья производная положения по времени)

Сложение этих двух расширений дает

$\mathbf {x} (t+\Delta t)=2\mathbf {x} (t)-\mathbf {x} (t-\Delta t)+\mathbf {a} (t)\Delta t^{2}+{\mathcal {O}}\left(\Delta t^{4}\right).$

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

Следует учитывать тот факт, что ускорение здесь вычисляется на основе точного решения: $\mathbf {a} (t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}$, тогда как на итерации он вычисляется в центральной точке итерации, $\mathbf {a} _{n}=\mathbf {A} (\mathbf {x} _{n})$. При вычислении глобальной ошибки, то есть расстояния между точным решением и последовательностью аппроксимации, эти два члена не сокращаются точно, что влияет на порядок глобальной ошибки.

Простой пример

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

Рассмотрим линейное дифференциальное уравнение ${\ddot {x}}(t)=w^{2}x(t)$ с постоянной $w$. Его точные базисные решения: $e^{wt}$ и $e^{-wt}$.

Метод Штёрмера, примененный к этому дифференциальному уравнению, приводит к линейному рекуррентному соотношению

$x_{n+1}-2x_{n}+x_{n-1}=h^{2}w^{2}x_{n},$

или

$x_{n+1}-2\left(1+{\tfrac {1}{2}}(wh)^{2}\right)x_{n}+x_{n-1}=0.$

Ее можно решить, найдя корни ее характеристического многочлена. $q^{2}-2\left(1+{\tfrac {1}{2}}(wh)^{2}\right)q+1=0$. Это

$q_{\pm }=1+{\tfrac {1}{2}}(wh)^{2}\pm wh{\sqrt {1+{\tfrac {1}{4}}(wh)^{2}}}.$

Базисные решения линейной рекуррентности: $x_{n}=q_{+}^{n}$ и $x_{n}=q_{-}^{n}$. Чтобы сравнить их с точными решениями, вычисляются разложения Тейлора:

${\begin{aligned}q_{+}&=1+{\tfrac {1}{2}}(wh)^{2}+wh\left(1+{\tfrac {1}{8}}(wh)^{2}-{\tfrac {3}{128}}(wh)^{4}+{\mathcal {O}}\left(h^{6}\right)\right)\\&=1+(wh)+{\tfrac {1}{2}}(wh)^{2}+{\tfrac {1}{8}}(wh)^{3}-{\tfrac {3}{128}}(wh)^{5}+{\mathcal {O}}\left(h^{7}\right).\end{aligned}}$

Частное этого ряда с показателем экспоненты $e^{wh}$ начинается с $1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)$, так

${\begin{aligned}q_{+}&=\left(1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)\right)e^{wh}\\&=e^{-{\frac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)}\,e^{wh}.\end{aligned}}$

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

${\begin{aligned}x_{n}=q_{+}^{n}&=e^{-{\frac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}\left(h^{4}\right)}\,e^{wt_{n}}\\&=e^{wt_{n}}\left(1-{\tfrac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}(h^{4})\right)\\&=e^{wt_{n}}+{\mathcal {O}}\left(h^{2}t_{n}e^{wt_{n}}\right).\end{aligned}}$

То есть, хотя локальная ошибка дискретизации имеет порядок 4, из-за второго порядка дифференциального уравнения глобальная ошибка имеет порядок 2 с константой, которая растет экспоненциально во времени.

Начинаем итерацию

Обратите внимание, что в начале итерации Verlet на шаге $n=1$, время $t=t_{1}=\Delta t$, вычисления $\mathbf {x} _{2}$, уже нужен вектор положения $\mathbf {x} _{1}$ во время $t=t_{1}$. На первый взгляд, это может вызвать проблемы, поскольку начальные условия известны только в начальный момент времени. $t_{0}=0$. Однако от этого ускорение $\mathbf {a} _{0}=\mathbf {A} (\mathbf {x} _{0})$ известно, и подходящую аппроксимацию положения на первом временном шаге можно получить с помощью полинома Тейлора второй степени:

$\mathbf {x} _{1}=\mathbf {x} _{0}+\mathbf {v} _{0}\Delta t+{\tfrac {1}{2}}\mathbf {a} _{0}\Delta t^{2}\approx \mathbf {x} (\Delta t)+{\mathcal {O}}\left(\Delta t^{3}\right).$

Тогда ошибка на первом временном шаге имеет порядок. ${\mathcal {O}}\left(\Delta t^{3}\right)$. Это не считается проблемой, поскольку при моделировании с большим количеством временных шагов ошибка на первом временном шаге составляет лишь пренебрежимо малую часть общей ошибки, которая во времени $t_{n}$ имеет порядок ${\mathcal {O}}\left(e^{Lt_{n}}\Delta t^{2}\right)$, как для расстояния векторов положения $\mathbf {x} _{n}$ к $\mathbf {x} (t_{n})$ что касается расстояния разделенных разностей ${\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}$ к ${\tfrac {\mathbf {x} (t_{n+1})-\mathbf {x} (t_{n})}{\Delta t}}$. Более того, чтобы получить эту глобальную ошибку второго порядка, первоначальная ошибка должна быть как минимум третьего порядка.

Непостоянные разницы во времени

Недостатком метода Штёрмера–Верле является то, что если шаг по времени ( $\Delta t$) изменяется, метод не аппроксимирует решение дифференциального уравнения. Это можно исправить по формуле

$\mathbf {x} _{i+1}=\mathbf {x} _{i}+\left(\mathbf {x} _{i}-\mathbf {x} _{i-1}\right){\frac {\Delta t_{i}}{\Delta t_{i-1}}}+\mathbf {a} _{i}\Delta t_{i}^{2}.$

Более точный вывод использует ряд Тейлора (до второго порядка) при $t_{i}$ на раз $t_{i+1}=t_{i}+\Delta t_{i}$ и $t_{i-1}=t_{i}-\Delta t_{i-1}$ получить после устранения $\mathbf {v} _{i}$

${\frac {\mathbf {x} _{i+1}-\mathbf {x} _{i}}{\Delta t_{i}}}+{\frac {\mathbf {x} _{i-1}-\mathbf {x} _{i}}{\Delta t_{i-1}}}=\mathbf {a} _{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}},$

так что итерационная формула принимает вид

$\mathbf {x} _{i+1}=\mathbf {x} _{i}+(\mathbf {x} _{i}-\mathbf {x} _{i-1}){\frac {\Delta t_{i}}{\Delta t_{i-1}}}+\mathbf {a} _{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}}\,\Delta t_{i}.$

метод Штёрмера –

Скорости не указаны явно в основном уравнении Штёрмера, но часто они необходимы для расчета определенных физических величин, таких как кинетическая энергия. Это может создать технические проблемы при моделировании молекулярной динамики , поскольку кинетическая энергия и мгновенные температуры во времени $t$ не может быть рассчитано для системы, пока позиции не известны во времени $t+\Delta t$. Этот недостаток можно устранить либо с помощью алгоритма Верле скорости , либо путем оценки скорости с использованием условий положения и теоремы о среднем значении :

$\mathbf {v} (t)={\frac {\mathbf {x} (t+\Delta t)-\mathbf {x} (t-\Delta t)}{2\Delta t}}+{\mathcal {O}}\left(\Delta t^{2}\right).$

Обратите внимание, что этот член скорости на шаг отстает от члена положения, поскольку он относится к скорости во времени. $t$, нет $t+\Delta t$, это означает, что $\mathbf {v} _{n}={\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n-1}}{2\Delta t}}$ является приближением второго порядка $\mathbf {v} (t_{n})$. Используя тот же аргумент, но сократив вдвое шаг по времени, $\mathbf {v} _{n+{\frac {1}{2}}}={\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}$ является приближением второго порядка $\mathbf {v} \left(t_{n+{\frac {1}{2}}}\right)$, с $t_{n+{\frac {1}{2}}}=t_{n}+{\tfrac {1}{2}}\Delta t$.

Можно сократить интервал, чтобы аппроксимировать скорость во времени $t+\Delta t$ ценой точности:

$\mathbf {v} (t+\Delta t)={\frac {\mathbf {x} (t+\Delta t)-\mathbf {x} (t)}{\Delta t}}+{\mathcal {O}}(\Delta t).$

Скорость Верле

Родственным и более часто используемым алгоритмом является алгоритм Верле скорости. аналогичен методу чехарды , за исключением того, что скорость и положение вычисляются при одном и том же значении переменной времени (чехарда этого не делает, как следует из названия). Здесь используется аналогичный подход, но явно учитывается скорость, решая проблему первого временного шага в базовом алгоритме Верле:

${\begin{aligned}\mathbf {x} (t+\Delta t)&=\mathbf {x} (t)+\mathbf {v} (t)\,\Delta t+{\tfrac {1}{2}}\,\mathbf {a} (t)\Delta t^{2},\\[6pt]\mathbf {v} (t+\Delta t)&=\mathbf {v} (t)+{\frac {\mathbf {a} (t)+\mathbf {a} (t+\Delta t)}{2}}\Delta t.\end{aligned}}$

Можно показать, что ошибка в скорости Верле того же порядка, что и в базовом Верле. Обратите внимание, что алгоритм скорости не обязательно требует больше памяти, поскольку в базовом Верле мы отслеживаем два вектора положения, а в Верле скорости мы отслеживаем один вектор положения и один вектор скорости. Стандартная схема реализации этого алгоритма:

1. Рассчитать $\mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)=\mathbf {v} (t)+{\tfrac {1}{2}}\,\mathbf {a} (t)\,\Delta t$.
2. Рассчитать $\mathbf {x} (t+\Delta t)=\mathbf {x} (t)+\mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)\,\Delta t$.
3. Вывести $\mathbf {a} (t+\Delta t)$ из потенциала взаимодействия, используя $\mathbf {x} (t+\Delta t)$.
4. Рассчитать $\mathbf {v} (t+\Delta t)=\mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)+{\tfrac {1}{2}}\,\mathbf {a} (t+\Delta t)\Delta t$.

Этот алгоритм также работает с переменными шагами по времени и идентичен форме интеграции метода чехарды «удар-дрейф-удар» .

Устранив полушаговую скорость, этот алгоритм можно сократить до

1. Рассчитать $\mathbf {x} (t+\Delta t)=\mathbf {x} (t)+\mathbf {v} (t)\,\Delta t+{\tfrac {1}{2}}\,\mathbf {a} (t)\,\Delta t^{2}$.
2. Вывести $\mathbf {a} (t+\Delta t)$ из потенциала взаимодействия, используя $\mathbf {x} (t+\Delta t)$.
3. Рассчитать $\mathbf {v} (t+\Delta t)=\mathbf {v} (t)+{\tfrac {1}{2}}\,{\bigl (}\mathbf {a} (t)+\mathbf {a} (t+\Delta t){\bigr )}\Delta t$.

Однако обратите внимание, что этот алгоритм предполагает, что ускорение $\mathbf {a} (t+\Delta t)$ зависит только от позиции $\mathbf {x} (t+\Delta t)$ и не зависит от скорости $\mathbf {v} (t+\Delta t)$.

Можно отметить, что долгосрочные результаты скорости Верле и, аналогично, чехарды на порядок лучше, чем полунеявный метод Эйлера . Алгоритмы практически идентичны с точностью до смещения скорости на полшага по времени. Это можно доказать, повернув приведенный выше контур, чтобы начать с шага 3, а затем заметив, что член ускорения на шаге 1 можно исключить, объединив шаги 2 и 4. Единственная разница состоит в том, что скорость в средней точке скорости Верле считается конечной скоростью. в полунеявном методе Эйлера.

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

Скоростной метод Верле является частным случаем бета-метода Ньюмарка с $\beta =0$ и $\gamma ={\tfrac {1}{2}}$.

представление

Поскольку скорость Verlet является обычно полезным алгоритмом в 3D-приложениях, общее решение, написанное на C++, может выглядеть, как показано ниже. Упрощенная сила сопротивления используется для демонстрации изменения ускорения, однако она необходима только в том случае, если ускорение непостоянно.

struct Body
{
    Vec3d pos { 0.0, 0.0, 0.0 };
    Vec3d vel { 2.0, 0.0, 0.0 }; // 2 m/s along x-axis
    Vec3d acc { 0.0, 0.0, 0.0 }; // no acceleration at first
    double mass = 1.0; // 1kg
    double drag = 0.1; // rho*C*Area – simplified drag for this example

    /**
     * Update pos and vel using "Velocity Verlet" integration
     * @param dt DeltaTime / time step [eg: 0.01]
     */
    void update(double dt)
    {
        Vec3d new_pos = pos + vel*dt + acc*(dt*dt*0.5);
        Vec3d new_acc = apply_forces(); // only needed if acceleration is not constant
        Vec3d new_vel = vel + (acc+new_acc)*(dt*0.5);
        pos = new_pos;
        vel = new_vel;
        acc = new_acc;
    }

    Vec3d apply_forces() const
    {
        Vec3d grav_acc = Vec3d{0.0, 0.0, -9.81 }; // 9.81 m/s² down in the z-axis
        Vec3d drag_force = 0.5 * drag * (vel * vel); // D = 0.5 * (rho * C * Area * vel^2)
        Vec3d drag_acc = drag_force / mass; // a = F/m
        return grav_acc - drag_acc;
    }
};

Условия ошибки

Основная статья: Ошибка усечения (численное интегрирование)

Глобальная ошибка усечения метода Верле равна ${\mathcal {O}}\left(\Delta t^{2}\right)$, как для положения, так и для скорости. Это контрастирует с тем фактом, что локальная ошибка в положении составляет всего лишь ${\mathcal {O}}\left(\Delta t^{4}\right)$ как описано выше. Разница связана с накоплением локальной ошибки усечения на всех итерациях.

Глобальную ошибку можно получить, приняв во внимание следующее:

$\operatorname {error} {\bigl (}x(t_{0}+\Delta t){\bigr )}={\mathcal {O}}\left(\Delta t^{4}\right)$

и

$x(t_{0}+2\Delta t)=2x(t_{0}+\Delta t)-x(t_{0})+\Delta t^{2}{\ddot {x}}(t_{0}+\Delta t)+{\mathcal {O}}\left(\Delta t^{4}\right).$

Поэтому

$\operatorname {error} {\bigl (}x(t_{0}+2\Delta t){\bigr )}=2\cdot \operatorname {error} {\bigl (}x(t_{0}+\Delta t){\bigr )}+{\mathcal {O}}\left(\Delta t^{4}\right)=3\,{\mathcal {O}}\left(\Delta t^{4}\right).$

Сходным образом:

${\begin{aligned}\operatorname {error} {\bigl (}x(t_{0}+3\Delta t){\bigl )}&=6\,{\mathcal {O}}\left(\Delta t^{4}\right),\\[6px]\operatorname {error} {\bigl (}x(t_{0}+4\Delta t){\bigl )}&=10\,{\mathcal {O}}\left(\Delta t^{4}\right),\\[6px]\operatorname {error} {\bigl (}x(t_{0}+5\Delta t){\bigl )}&=15\,{\mathcal {O}}\left(\Delta t^{4}\right),\end{aligned}}$

которое можно обобщить до (это можно показать по индукции, но здесь оно приводится без доказательства):

$\operatorname {error} {\bigl (}x(t_{0}+n\Delta t){\bigr )}={\frac {n(n+1)}{2}}\,{\mathcal {O}}\left(\Delta t^{4}\right).$

Если мы рассмотрим глобальную ошибку в положении между $x(t)$ и $x(t+T)$, где $T=n\Delta t$, ясно, что

$\operatorname {error} {\bigl (}x(t_{0}+T){\bigr )}=\left({\frac {T^{2}}{2\Delta t^{2}}}+{\frac {T}{2\Delta t}}\right){\mathcal {O}}\left(\Delta t^{4}\right),$

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

$\operatorname {error} {\bigr (}x(t_{0}+T){\bigl )}={\mathcal {O}}\left(\Delta t^{2}\right).$

Поскольку скорость определяется некумулятивным способом по позициям интегратора Верле, глобальная ошибка скорости также равна ${\mathcal {O}}\left(\Delta t^{2}\right)$.

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

Ограничения

Основная статья: Алгоритм ограничений

Системы множества частиц с ограничениями проще решать с помощью интегрирования Верле, чем с помощью методов Эйлера. Ограничениями между точками могут быть, например, потенциалы, удерживающие их на определенном расстоянии, или силы притяжения. Их можно смоделировать как пружины, соединяющие частицы. Используя пружины бесконечной жесткости, модель можно затем решить с помощью алгоритма Верле.

В одном измерении взаимосвязь между неограниченными позициями ${\tilde {x}}_{i}^{(t)}$ и реальные позиции $x_{i}^{(t)}$ очков $i$ во время $t$, учитывая желаемое расстояние ограничения $r$, можно найти по алгоритму

${\begin{aligned}d_{1}&=x_{2}^{(t)}-x_{1}^{(t)},\\[6px]d_{2}&=\|d_{1}\|,\\[6px]d_{3}&={\frac {d_{2}-r}{d_{2}}},\\[6px]x_{1}^{(t+\Delta t)}&={\tilde {x}}_{1}^{(t+\Delta t)}+{\tfrac {1}{2}}d_{1}d_{3},\\[6px]x_{2}^{(t+\Delta t)}&={\tilde {x}}_{2}^{(t+\Delta t)}-{\tfrac {1}{2}}d_{1}d_{3}.\end{aligned}}$

Интеграция Верле полезна, поскольку она напрямую связывает силу с положением, а не решает проблему с использованием скоростей.

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

При локальной аппроксимации ограничений до первого порядка это то же самое, что и метод Гаусса – Зейделя . для небольших матриц Известно, что LU-разложение происходит быстрее. Большие системы можно разделить на кластеры (например, каждая тряпичная кукла = кластер). Внутри кластеров используется метод LU, между кластерами — метод Гаусса–Зейделя . Матричный код можно использовать повторно: зависимость сил от положений можно локально аппроксимировать до первого порядка, а интегрирование Верле можно сделать более неявным.

Сложное программное обеспечение, такое как SuperLU существует для решения сложных задач с использованием разреженных матриц. Конкретные методы, такие как использование (кластеров) матриц, могут использоваться для решения конкретной проблемы, например, проблемы распространения силы через лист ткани без образования звуковой волны .

Другой способ решения голономных ограничений — использование алгоритмов ограничений .

Реакции на столкновение

Одним из способов реагирования на столкновения является использование системы штрафов, которая по сути применяет заданную силу к точке при контакте. Проблема в том, что очень сложно выбрать прикладываемую силу. Используйте слишком большую силу, и объекты станут нестабильными, слишком слабыми, и объекты проникнут друг в друга. Другой способ — использовать реакции столкновения проекций, которые берут точку нарушения и пытаются переместить ее на кратчайшее возможное расстояние, чтобы удалить ее от другого объекта.

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

Два простейших метода определения новой скорости — это совершенно упругие и неупругие столкновения . Немного более сложная стратегия, предлагающая больший контроль, предполагает использование коэффициента реституции .

См. также

• Условие Куранта – Фридрихса – Леви.
• Энергетический дрейф
• Симплектический интегратор
• Чехарда интеграции
• Алгоритм Бимана

Внешние ссылки

• Демонстрация интеграции Verlet и код в виде Java-апплета
• Продвинутая физика персонажей Томаса Якобсена
• Теория молекулярно-динамического моделирования – внизу страницы
• Интеграция Verlet реализована на современном JavaScript – внизу страницы.