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

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

• 1 Базовый Störmer – Verlet
  • 1.1 Уравнения движения
  • 1.2 Интегрирование Верле (без скоростей)
  • 1.3 Ошибка дискретизации
  • 1.4 Простой пример
  • 1.5 Начало итерации
  • 1.6 Непостоянная разница во времени
  • 1.7 Вычисление скоростей - метод Штёрмера – Верле
• 2 Velocity Verlet
  • 2.1 Алгоритмическое представление
• 3 Термины ошибок
• 4 Ограничения
• 5 Реакции коллизий
• 6 См. Также
• 7 Литература
• 8 Внешние ссылки

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

• установить $x\; \to \; 1\; =\; x\; \to \; 0\; +\; v\; \to \; 0\; \Delta \; T\; +\; 1\; 2\; A\; \to \; (Икс\; \to \; 0)\; \Delta \; T\; 2\; \{\backslash \; Displaystyle\; \{\backslash \; vec\; \{x\}\}\; \_\; \{1\}\; =\; \{\backslash \; vec\; \{x\}\}\; \_\; \{0\}\; +\; \{\backslash \; vec\; \{v\}\}\; \_\; \{0\}\; \backslash ,\; \backslash \; Delta\; t\; +\; \{\backslash \; frac\; \{1\}\; \{2\}\}\; \{\backslash \; vec\; \{A\}\}\; (\{\backslash \; vec\; \{x\}\}\; \_\; \{0\})\; \backslash ,\; \backslash \; Delta\; t\; ^\; \{2\}\}$,
• для n = 1, 2,... итерация

  $x\; \to \; n\; +\; 1\; =\; 2\; x\; \to \; n\; -\; x\; \to \; n\; -\; 1\; +\; A\; \to \; (x\; \to \; n)\; \Delta \; t\; 2.\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n\; +\; 1\}\; =\; 2\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n\}\; -\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n-1\}\; +\; \{\backslash \; vec\; \{A\}\; \}\; (\{\backslash \; vec\; \{x\}\}\; \_\; \{n\})\; \backslash ,\; \backslash \; Delta\; t\; ^\; \{2\}.\}$

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

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

$M\; Икс\; \to \; ¨\; (T)\; знак\; равно\; F\; (Икс\; \to \; (T))\; =\; -\; \nabla \; V\; (x\; \to \; (T)),\; \{\backslash \; Displaystyle\; M\; \{\backslash \; ddot\; \{\backslash \; vec\; \{x\}\}\}\; (т)\; =\; F\; \{\backslash \; big\; (\}\; \{\backslash \; vec\; \{x\}\}\; (t)\; \{\backslash \; big)\}\; =\; -\; \backslash \; nabla\; V\; \{\backslash \; big\; (\}\; \{\backslash \; vec\; \{x\}\}\; (t)\; \{\backslash \; big)\},\}$

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

$mkx\; \to \; ¨\; К\; (T)\; =\; F\; К\; (Икс\; \to \; (T))\; =\; -\; \nabla \; Икс\; \to \; К\; V\; (X\; \to \; (T)),\; \{\backslash \; Displaystyle\; M\_\; \{k\}\; \{\backslash \; ddot\; \{\backslash \; vec\; \{x\}\; \}\}\; \_\; \{k\}\; (t)\; =\; F\_\; \{k\}\; \{\backslash \; big\; (\}\; \{\backslash \; vec\; \{x\}\}\; (t)\; \{\backslash \; big)\}\; =\; -\; \backslash \; nabla\; \_\; \{\{\backslash \; vec\; \{x\}\}\; \_\; \{k\; \}\}\; V\; \{\backslash \; big\; (\}\; \{\backslash \; vec\; \{x\}\}\; (t)\; \{\backslash \; big)\},\}$

где

t - время,

$x\; \to \; (t)\; =\; (x\; \to \; 1\; (т),\dots ,\; Икс\; \to \; N\; (T))\; \{\backslash \; Displaystyle\; \{\backslash \; vec\; \{x\}\}\; (т)\; =\; \{\backslash \; big\; (\}\; \{\backslash \; vec\; \{x\}\}\; \_\; \{1\}\; (т),\; \backslash \; ldots,\; \{\backslash \; vec\; \{x\}\}\; \_\; \{N\}\; (t)\; \{\backslash \; big)\}\}$- это ансамбль вектора положения N объектов,

V - скалярная потенциальная функция,

F - отрицательный градиент потенциала, задающий ансамбль сил, действующих на частицу. s,

M - матрица масс, обычно диагональная с блоками с массой $mk\; \{\backslash \; displaystyle\; m\_\; \{k\}\}$для каждой частицы.

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

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

$x\; \to \; ¨\; (t)\; =\; A\; (x\; \to \; (t))\; \{\backslash \; displaystyle\; \{\backslash \; ddot\; \{\backslash \; vec\; \{x\}\}\}\; (t)\; =\; A\; \{\backslash \; big\; (\}\; \{\backslash \; vec\; \{x\}\}\; (t)\; \{\backslash \; big)\}\}$

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

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

Для дискретизации и численного решения этой задачи начального значения временной шаг $\Delta \; t>0\; \{\backslash \; displaystyle\; \backslash \; Delta\; t>0\}$выбирается, и рассматривается последовательность точек выборки $tn\; =\; n\; \Delta \; t\; \{\backslash \; displaystyle\; t\_\; \{n\}\; =\; n\; \backslash ,\; \backslash \; Delta\; t\}$. Задача состоит в том, чтобы построить последовательность точки $x\; \to \; n\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n\}\}$, которые близко следуют за точками $x\; \to \; (tn)\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\; \}\; (t\_\; \{n\})\}$на траектории точного решения.

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

$\Delta \; 2\; x\; \to \; n\; \Delta \; t\; 2\; =\; x\; \to \; n\; +\; 1\; -\; x\; \to \; n\; \Delta \; t\; -\; x\; \to \; n\; -\; x\; \to \; n\; -\; 1\; \Delta \; t\; \Delta \; t\; =\; x\; \to \; n\; +\; 1\; -\; 2\; x\; \to \; n\; +\; x\; \to \; n\; -\; 1\; \Delta \; t\; 2\; =\; a\; \to \; n\; =\; A\; \to \; (x\; \to \; n).\; \{\backslash \; displaystyle\; \{\backslash \; frac\; \{\backslash \; Delta\; ^\; \{2\}\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n\}\}\; \{\backslash \; Delta\; t\; ^\; \{2\}\}\}\; =\; \{\backslash \; frac\; \{\{\backslash \; frac\; \{\{\backslash \; vec\; \{x\})\; \}\; \_\; \{n\; +\; 1\}\; -\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n\}\}\; \{\backslash \; Delta\; t\}\}\; -\; \{\backslash \; frac\; \{\{\backslash \; vec\; \{x\}\}\; \_\; \{n\}\; -\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n-1\}\}\; \{\backslash \; Delta\; t\}\}\}\; \{\backslash \; Delta\; t\}\}\; =\; \{\backslash \; frac\; \{\{\backslash \; vec\; \{x\}\}\; \_\; \{n\; +\; 1\}\; -2\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n\; \}\; +\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n-1\}\}\; \{\backslash \; Delta\; t\; ^\; \{2\}\}\}\; =\; \{\backslash \; vec\; \{a\}\}\; \_\; \{n\}\; =\; \{\backslash \; vec\; \{A\}\}\; (\{\backslash \; vec\; \{x\}\}\; \_\; \{n\}).\}$

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

$x\; \to \; n\; +\; 1\; =\; 2\; x\; \to \; n\; -\; x\; \to \; n\; -\; 1\; +\; a\; \to \; n\; \Delta \; t\; 2,\; a\; \to \; n\; =\; A\; \to \; (x\; \to \; n).\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n\; +\; 1\}\; =\; 2\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n\}\; -\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n-1\}\; +\; \{\backslash \; vec\; \{a\}\; \}\; \_\; \{n\}\; \backslash ,\; \backslash \; Delta\; t\; ^\; \{2\},\; \backslash \; qquad\; \{\backslash \; vec\; \{a\}\}\; \_\; \{n\}\; =\; \{\backslash \; vec\; \{A\}\}\; (\{\backslash \; vec\; \{x\}\}\; \_\; \{n\}).\; \}$

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

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

$x\; \to \; (t\; +\; \Delta \; t)\; =\; x\; \to \; (t)\; +\; v\; \to \; (t)\; \Delta \; t\; +\; a\; \to \; (t)\; \Delta \; t\; 2\; 2\; +\; b\; \to \; (t)\; \Delta \; t\; 3\; 6\; +\; O\; (\Delta \; t\; 4)\; x\; \to \; (t\; -\; \Delta \; t)\; =\; x\; \to \; (t)\; -\; v\; \to \; (t)\; \Delta \; t\; +\; a\; \to \; (\; t)\; \Delta \; t\; 2\; 2\; -\; b\; \to \; (t)\; \Delta \; t\; 3\; 6\; +\; O\; (\Delta \; t\; 4),\; \{\backslash \; displaystyle\; \{\backslash \; begin\; \{align\}\; \{\backslash \; vec\; \{x\}\}\; (t\; +\; \backslash \; Delta\; t)\; =\; \{\backslash \; vec\; \{x\}\}\; (t)\; +\; \{\backslash \; vec\; \{v\}\}\; (t)\; \backslash \; Delta\; t\; +\; \{\backslash \; frac\; \{\{\backslash \; vec\; \{a\}\}\; (t)\; \backslash \; Delta\; t\; ^\; \{2\}\}\; \{2\}\}\; +\; \{\; \backslash \; frac\; \{\{\backslash \; vec\; \{b\}\}\; (t)\; \backslash \; Delta\; t\; ^\; \{3\}\}\; \{6\}\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (\backslash \; Delta\; t\; ^\; \{4\})\; \backslash \backslash \; \{\backslash \; vec\; \{x\}\; \}\; (t-\; \backslash \; Delta\; t)\; =\; \{\backslash \; vec\; \{x\}\}\; (t)\; -\; \{\backslash \; vec\; \{v\}\}\; (t)\; \backslash \; Delta\; t\; +\; \{\backslash \; frac\; \{\{\backslash \; vec\; \{a\}\}\; (t)\; \backslash \; Delta\; t\; ^\; \{2\}\}\; \{2\}\}\; -\; \{\backslash \; frac\; \{\{\backslash \; vec\; \{b\}\}\; (t)\; \backslash \; Delta\; t\; ^\; \{3\}\}\; \{6\}\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (\backslash \; De\; lta\; t\; ^\; \{4\}),\; \backslash \; end\; \{align\}\}\}$

где $x\; \to \; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\}$- позиция, $v\; \to \; =\; x\; \to \; \dot{}\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\}\}\; =\; \{\backslash \; dot\; \{\backslash \; vec\; \{x\}\}\}\}$скорость, $a\; \to \; =\; x\; \to \; ¨\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{a\}\}\; =\; \{\backslash \; ddot\; \{\backslash \; vec\; \{x\}\}\}\}$ускорение и $b\; \to \; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{b\}\}\}$рывок (третья производная положения по времени).

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

$x\; \to \; (t\; +\; \Delta \; t)\; =\; 2\; x\; \to \; (t)\; -\; x\; \to \; (t\; -\; \Delta \; t)\; +\; a\; \to \; (t)\; \Delta \; t\; 2\; +\; O\; (\; \Delta \; t\; 4).\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; (t\; +\; \backslash \; Delta\; t)\; =\; 2\; \{\backslash \; vec\; \{x\}\}\; (t)\; -\; \{\backslash \; vec\; \{x\}\}\; (t-\; \backslash \; Delta\; t)\; +\; \{\backslash \; vec\; \{a\}\; \}\; (t)\; \backslash \; Delta\; t\; ^\; \{2\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (\backslash \; Delta\; t\; ^\; \{4\}).\}$

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

Следует с осторожностью относиться к тому факту, что ускорение здесь вычисляется из точного решения, $a\; \to \; (t)\; =\; A\; (x\; \to \; (t))\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{a\}\; \}\; (t)\; =\; A\; \{\backslash \; big\; (\}\; \{\backslash \; vec\; \{x\}\}\; (t)\; \{\backslash \; big)\}\}$, тогда как в итерации он вычисляется в центральной точке итерации, $a\; \to \; N\; =\; A\; (x\; \to \; n)\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{a\}\}\; \_\; \{n\}\; =\; A\; (\{\backslash \; vec\; \{x\}\}\; \_\; \{n\})\}$. При вычислении глобальной ошибки, то есть расстояния между точным решением и аппроксимационной последовательностью, эти два члена не компенсируются точно, что влияет на порядок глобальной ошибки.

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

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

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

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

$xn\; +\; 1-2\; xn\; +\; xn\; -\; 1\; =\; h\; 2\; w\; 2\; xn,\; \{\backslash \; displaystyle\; x\_\; \{n\; +\; 1\}\; -2x\_\; \{n\}\; +\; x\_\; \{n-\; 1\}\; =\; h\; ^\; \{2\}\; w\; ^\; \{2\}\; x\_\; \{n\},\}$

или

$xn\; +\; 1-2\; (1\; +\; 1\; 2\; (wh)\; 2)\; xn\; +\; xn\; -\; 1\; =\; 0.\; \{\; \backslash \; Displaystyle\; x\_\; \{n\; +\; 1\}\; -2\; \backslash \; left\; (1\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; (wh)\; ^\; \{2\}\; \backslash \; right)\; x\_\; \{n\}\; +\; x\_\; \{n-1\}\; =\; 0.\}$

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

$q\; \pm \; =\; 1\; +\; 1\; 2\; (w\; h)\; 2\; \pm \; w\; h\; 1\; +\; 1\; 4\; (w\; h)\; 2.\; \{\backslash \; displaystyle\; q\; \_\; \{\backslash \; pm\}\; =\; 1\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; (wh)\; ^\; \{2\}\; \backslash \; pm\; wh\; \{\backslash \; sqrt\; \{1\; +\; \{\backslash \; tfrac\; \{1\}\; \{4\}\}\; (wh)\; ^\; \{2\}\}\}.\}$

Базовые решения линейной повторяемости: $xn\; =\; q\; +\; n\; \{\backslash \; displaystyle\; x\_\; \{n\}\; =\; q\; \_\; \{+\}\; ^\; \{n\}\}$и $xn\; =\; q\; -\; n\; \{\backslash \; displaystyle\; x\_\; \{n\}\; =\; q\; \_\; \{-\}\; ^\; \{n\}\}$. Чтобы сравнить их с точными решениями, вычисляются разложения Тейлора:

$q\; +\; =\; 1\; +\; 1\; 2\; (wh)\; 2\; +\; wh\; (1\; +\; 1\; 8\; (wh)\; 2\; -\; 3\; 128\; (wh)\; 4\; +\; O\; (h\; 6))\; =\; 1\; +\; (белый)\; +\; 1\; 2\; (белый)\; 2\; +\; 1\; 8\; (белый)\; 3\; -\; 3\; 128\; (белый)\; 5\; +\; O\; (h\; 7).\; \{\backslash \; displaystyle\; \{\backslash \; begin\; \{align\}\; q\; \_\; \{+\}\; =\; 1\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; (wh)\; ^\; \{2\}\; +\; wh\; \backslash \; left\; (1\; +\; \{\backslash \; tfrac\; \{1\}\; \{8\})\; \}\; (wh)\; ^\; \{2\}\; -\; \{\backslash \; tfrac\; \{3\}\; \{128\}\}\; (wh)\; ^\; \{4\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (h\; ^\; \{6\})\; \backslash \; right)\; \backslash \backslash \; =\; 1+\; (wh)\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; (wh)\; ^\; \{2\}\; +\; \{\backslash \; tfrac\; \{1\}\; \{8\}\}\; (wh)\; ^\; \{3\}\; -\; \{\backslash \; tfrac\; \{3\}\; \{128\}\}\; (wh)\; ^\; \{5\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (h\; ^\; \{7\}).\; \backslash \; end\; \{align\}\}\}$

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

$q\; +\; =\; (1\; -\; 1\; 24\; (wh)\; 3\; +\; O\; (h\; 5))\; ewh\; =\; e\; -\; 1\; 24\; (wh)\; 3\; +\; O\; (h\; 5)\; ewh.\; \{\backslash \; displaystyle\; \{\backslash \; begin\; \{align\}\; q\; \_\; \{+\}\; =\; \backslash \; left\; (1\; -\; \{\backslash \; tfrac\; \{1\}\; \{24\}\}\; (wh)\; ^\; \{3\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (h\; ^\; \{5\; \})\; \backslash \; right)\; e\; ^\; \{wh\}\; \backslash \backslash \; =\; e\; ^\; \{-\; \{\backslash \; frac\; \{1\}\; \{24\}\}\; (wh)\; ^\; \{3\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (h\; ^\; \{5\})\; \}\; \backslash ,\; e\; ^\; \{wh\}.\; \backslash \; end\; \{align\}\}\}$

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

$xn\; =\; q\; +\; n\; =\; e\; -\; 1\; 24\; (wh)\; 2\; wtn\; +\; O\; (h\; 4)\; ewtn\; =\; ewtn\; (1\; -\; 1\; 24\; (wh)\; 2\; wtn\; +\; O\; (h\; 4))\; =\; ewtn\; +\; O\; (h\; 2\; tnewtn).\; \{\backslash \; displaystyle\; \{\backslash \; begin\; \{align\}\; x\_\; \{n\}\; =\; q\; \_\; \{+\}\; ^\; \{\backslash ;\; n\}\; =\; e\; ^\; \{-\; \{\backslash \; frac\; \{1\}\; \{24\}\}\; (wh)\; ^\; \{2\}\; \backslash ,\; wt\_\; \{n\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (h\; ^\; \{4\})\}\; \backslash ,\; e\; ^\; \{wt\_\; \{n\}\}\; \backslash \backslash \; =\; e\; ^\; \{wt\_\; \{n\}\}\; \backslash \; left\; (1\; -\; \{\backslash \; tfrac\; \{\; 1\}\; \{24\}\}\; (wh)\; ^\; \{2\}\; \backslash ,\; wt\_\; \{n\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (h\; ^\; \{4\})\; \backslash \; right)\; \backslash \backslash \; =\; e\; ^\; \{wt\_\; \{n\}\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (h\; ^\; \{2\}\; t\_\; \{n\}\; e\; ^\; \{wt\_\; \{n\}\}).\; \backslash \; end\; \{align\}\}\}$

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

Запуск итерации

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

$x\; \to \; 1\; =\; x\; \to \; 0\; +\; v\; \to \; 0\; \Delta \; t\; +\; 1\; 2\; a\; \to \; 0\; \Delta \; t\; 2\; \approx \; x\; \to \; (\Delta \; t)\; +\; O\; (\Delta \; t\; 3).\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; \_\; \{1\}\; =\; \{\backslash \; vec\; \{x\}\}\; \_\; \{0\}\; +\; \{\backslash \; vec\; \{v\}\}\; \_\; \{0\}\; \backslash \; Delta\; t\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\; \}\; \{\backslash \; vec\; \{a\}\}\; \_\; \{0\}\; \backslash \; Delta\; t\; ^\; \{2\}\; \backslash \; приблизительно\; \{\backslash \; vec\; \{x\}\}\; (\backslash \; Delta\; t)\; +\; \{\backslash \; mathcal\; \{O\}\}\; (\backslash \; Delta\; t\; ^\; \{3\}).\}$

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

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

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

$x\; \to \; i\; +\; 1\; =\; x\; \to \; i\; +\; (x\; \to \; i\; -\; x\; \to \; i\; -\; 1)\; (\Delta \; t\; i\; /\; \Delta \; t\; i\; -\; 1)\; +\; a\; \to \; i\; \Delta \; t\; i\; 2.\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; \_\; \{я\; +\; 1\}\; =\; \{\backslash \; vec\; \{x\}\}\; \_\; \{i\}\; +\; (\{\backslash \; vec\; \{x\}\}\; \_\; \{i\}\; -\; \{\backslash \; vec\; \{x\}\}\; \_\; \{i-1\})\; (\backslash \; Delta\; t\_\; \{i\}\; /\; \backslash \; Delta\; t\_\; \{i-1\})\; +\; \{\backslash \; vec\; \{a\}\}\; \_\; \{i\}\; \backslash \; Delta\; t\_\; \{i\}\; ^\; \{2\}.\}$

Более точный вывод использует ряд Тейлора (до второго порядка) в $ti\; \{\backslash \; displaystyle\; t\_\; \{i\}\}$для времен $ti\; +\; 1\; =\; ti\; +\; \Delta \; ti\; \{\backslash \; displaystyle\; t\_\; \{\; я\; +\; 1\}\; =\; t\_\; \{i\}\; +\; \backslash \; Delta\; t\_\; \{i\}\}$и $ti\; -\; 1\; =\; ti\; -\; \Delta \; ti\; -\; 1\; \{\backslash \; displaystyle\; t\_\; \{i-1\}\; =\; t\_\; \{i\; \}\; -\; \backslash \; Delta\; t\_\; \{i-1\}\}$для получения после исключения $v\; \to \; i\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\}\}\; \_\; \{i\}\}$

$x\; \to \; i\; +\; 1\; -\; Икс\; \to \; я\; \Delta \; ti\; +\; x\; \to \; я\; -\; 1\; -\; x\; \to \; я\; \Delta \; ti\; -\; 1\; =\; a\; \to \; я\; \Delta \; ti\; +\; \Delta \; ti\; -\; 1\; 2,\; \{\backslash \; displaystyle\; \{\backslash \; frac\; \{\{\backslash \; vec\; \{x\}\}\; \_\; \{\; i\; +\; 1\}\; -\; \{\backslash \; vec\; \{x\}\}\; \_\; \{i\}\}\; \{\backslash \; Delta\; t\_\; \{i\}\}\}\; +\; \{\backslash \; frac\; \{\{\backslash \; vec\; \{x\}\}\; \_\; \{i-1\}\; -\; \{\backslash \; vec\; \{x\; \}\}\; \_\; \{i\}\}\; \{\backslash \; Delta\; t\_\; \{i-1\}\}\}\; =\; \{\backslash \; vec\; \{a\}\}\; \_\; \{i\}\; \backslash ,\; \{\backslash \; frac\; \{\backslash \; Delta\; t\_\; \{i\}\; +\; \backslash \; Delta\; t\_\; \{i-1\; \}\}\; \{2\}\},\}$

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

$x\; \to \; i\; +\; 1\; =\; x\; \to \; i\; +\; (x\; \to \; i\; -\; x\; \to \; i\; -\; 1)\; \Delta \; ti\; \Delta \; ti\; -\; 1\; +\; a\; \to \; i\; \Delta \; ti\; +\; \Delta \; ti\; -\; 1\; 2\; \Delta \; ti.\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; \_\; \{я\; +\; 1\}\; =\; \{\backslash \; vec\; \{x\}\}\; \_\; \{i\}\; +\; (\{\backslash \; vec\; \{x\}\}\; \_\; \{i\}\; -\; \{\backslash \; vec\; \{x\}\}\; \_\; \{i-1\})\; \{\backslash \; frac\; \{\backslash \; Delta\; t\_\; \{i\}\}\; \{\backslash \; Delta\; t\_\; \{i-1\}\}\}\; +\; \{\backslash \; vec\; \{a\}\}\; \_\; \{i\}\; \backslash ,\; \{\backslash \; frac\; \{\backslash \; Delta\; t\_\; \{\; i\}\; +\; \backslash \; Delta\; t\_\; \{i-1\}\}\; \{2\}\}\; \backslash ,\; \backslash \; Delta\; t\_\; \{i\}.\}$

Вычисление скоростей - метод Штёрмера – Верле

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

$v\; \to \; (t)\; =\; x\; \to \; (t\; +\; \Delta \; t)\; -\; x\; \to \; (t\; -\; \Delta \; t)\; 2\; \Delta \; t\; +\; O\; (\Delta \; t\; 2).\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\}\}\; (t)\; =\; \{\backslash \; frac\; \{\{\backslash \; vec\; \{x\}\}\; (t\; +\; \backslash \; Delta\; t)\; -\; \{\backslash \; vec\; \{x\}\}\; (t-\; \backslash \; Delta\; t)\}\; \{2\; \backslash \; Delta\; t\}\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (\backslash \; Delta\; t\; ^\; \{2\}).\}$

Обратите внимание, что этот член скорости на шаг отстает от члена положения, поскольку он предназначен для скорости в момент времени $t\; \{\backslash \; displaystyle\; t\}$, а не $t\; +\; \Delta \; t\; \{\backslash \; displaystyle\; t\; +\; \backslash \; Delta\; t\}$, что означает, что $v\; \to \; n\; =\; x\; \to \; n\; +\; 1\; -\; Икс\; \to \; N\; -\; 1\; 2\; \Delta \; T\; \{\backslash \; Displaystyle\; \{\backslash \; vec\; \{v\}\}\; \_\; \{n\}\; =\; \{\backslash \; tfrac\; \{\{\backslash \; vec\; \{x\}\}\; \_\; \{n\; +\; 1\}\; -\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n-1\}\}\; \{2\; \backslash \; Delta\; t\}\}\}$- это приближение второго порядка к $v\; \to \; (tn)\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\}\}\; (t\_\; \{n\})\}$. Используя тот же аргумент, но с уменьшением вдвое временного шага, $v\; \to \; n\; +\; 1/2\; =\; x\; \to \; n\; +\; 1\; -\; x\; \to \; n\; \Delta \; t\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\}\}\; \_\; \{n\; +\; 1/2\; \}\; =\; \{\backslash \; tfrac\; \{\{\backslash \; vec\; \{x\}\}\; \_\; \{n\; +\; 1\}\; -\; \{\backslash \; vec\; \{x\}\}\; \_\; \{n\}\}\; \{\backslash \; Delta\; t\}\}\}$является второстепенным приближение к $v\; \to \; (tn\; +\; 1/2)\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\}\}\; (t\_\; \{n\; +\; 1/2\})\}$, с $tn\; +\; 1\; /\; 2\; =\; tn\; +\; 1\; 2\; \Delta \; t\; \{\backslash \; displaystyle\; t\_\; \{n\; +\; 1/2\}\; =\; t\_\; \{n\}\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; \backslash \; Delta\; t\}$.

Можно сократить интервал, чтобы приблизительно скорость во время $t\; +\; \Delta \; t\; \{\backslash \; displaystyle\; t\; +\; \backslash \; Delta\; t\}$за счет точности:

$v\; \to \; (t\; +\; \Delta \; t)\; =\; x\; \to \; (t\; +\; \Delta \; t)\; -\; x\; \to \; (t)\; \Delta \; t\; +\; O\; (\Delta \; t).\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\}\}\; (t\; +\; \backslash \; Delta\; t)\; =\; \{\backslash \; frac\; \{\{\backslash \; vec\; \{x\}\}\; (t\; +\; \backslash \; Delta\; t)\; -\; \{\backslash \; vec\; \{x\}\}\; (t)\}\; \{\backslash \; Delta\; t\; \}\}\; +\; \{\backslash \; mathcal\; \{O\}\}\; (\backslash \; Delta\; t).\}$

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

$x\; \to \; (t\; +\; \Delta \; t)\; =\; x\; \to \; (t)\; +\; v\; \to \; (t)\; \Delta \; t\; +\; 1\; 2\; a\; \to \; (t)\; \Delta \; t\; 2,\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; (t\; +\; \backslash \; Delta\; t)\; =\; \{\backslash \; vec\; \{x\}\}\; (t)\; +\; \{\backslash \; vec\; \{v\}\}\; (t)\; \backslash ,\; \backslash \; Delta\; t\; +\; \{\backslash \; frac\; \{1\}\; \{2\}\}\; \backslash ,\; \{\backslash \; vec\; \{a\}\}\; (t)\; \backslash \; Delta\; t\; ^\; \{2\},\}$

$v\; \to \; (t\; +\; \Delta \; t)\; =\; v\; \to \; (t)\; +\; a\; \to \; (t)\; +\; a\; \to \; (t\; +\; \Delta \; t)\; 2\; \Delta \; t.\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\}\}\; (t\; +\; \backslash \; Delta\; t)\; =\; \{\backslash \; vec\; \{v\}\}\; (t)\; +\; \{\backslash \; frac\; \{\{\backslash \; vec\; \{a\}\}\; (t)\; +\; \{\backslash \; vec\; \{a\}\}\; (t\; +\; \backslash \; Delta\; t)\}\; \{2\}\}\; \backslash \; Delta\; t.\}$

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

1. Вычислить $v\; \to \; (t\; +\; 1\; 2\; \Delta \; t)\; =\; v\; \to \; (t)\; +\; 1\; 2\; a\; \to \; (t)\; \Delta \; t\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\; \}\}\; \backslash \; left\; (t\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; \backslash ,\; \backslash \; Delta\; t\; \backslash \; right)\; =\; \{\backslash \; vec\; \{v\}\}\; (t)\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; \backslash ,\; \{\backslash \; vec\; \{a\}\}\; (t)\; \backslash ,\; \backslash \; Delta\; t\}$.
2. Вычислить $x\; \to \; (t\; +\; \Delta \; t)\; =\; x\; \to \; (t)\; +\; v\; \to \; (t\; +\; 1\; 2\; \Delta \; t)\; \Delta \; t\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; (t\; +\; \backslash \; Delta\; t)\; =\; \{\backslash \; vec\; \{x\}\}\; (t)\; +\; \{\backslash \; vec\; \{v\}\}\; \backslash \; left\; (t\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; \backslash ,\; \backslash \; Delta\; t\; \backslash \; right)\; \backslash ,\; \backslash \; Delta\; t\}$.
3. Вывести $a\; \to \; (t\; +\; \Delta \; t)\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{a\}\}\; (t\; +\; \backslash \; Delta\; t)\}$из потенциал взаимодействия с помощью $x\; \to \; (t\; +\; \Delta \; t)\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; (t\; +\; \backslash \; Delta\; t)\}$.
4. Вычислить $v\; \to \; (t\; +\; \Delta \; t)\; =\; v\; \to \; (T\; +\; 1\; 2\; \Delta \; t)\; +\; 1\; 2\; a\; \to \; (t\; +\; \Delta \; t)\; \Delta \; t\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\}\}\; (t\; +\; \backslash \; Delta\; t)\; =\; \{\backslash \; vec\; \{v\}\}\; \backslash \; left\; (t\; +\; \{\; \backslash \; tfrac\; \{1\}\; \{2\}\}\; \backslash ,\; \backslash \; Delta\; t\; \backslash \; right)\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; \backslash ,\; \{\backslash \; vec\; \{a\}\}\; (t\; +\; \backslash \; Delta\; t)\; \backslash \; Delta\; t\}$.

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

1. Вычислить $x\; \to \; (t\; +\; \Delta \; t)\; =\; x\; \to \; (t)\; +\; v\; \to \; (t)\; \Delta \; t\; +\; 1\; 2\; a\; \to \; (t)\; \Delta \; t\; 2\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; (t\; +\; \backslash \; Delta\; t)\; =\; \{\backslash \; vec\; \{x\}\}\; (t)\; +\; \{\backslash \; vec\; \{v\}\}\; (t)\; \backslash ,\; \backslash \; Delta\; t\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; \backslash ,\; \{\backslash \; vec\; \{a\}\}\; (t)\; \backslash ,\; \backslash \; Delta\; t\; ^\; \{2\}\}$.
2. Вывести $a\; \to \; (t\; +\; \Delta \; t)\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{a\}\; \}\; (t\; +\; \backslash \; Delta\; t)\}$из потенциала взаимодействия с использованием $x\; \to \; (t\; +\; \Delta \; t)\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{x\}\}\; (t\; +\; \backslash \; Delta\; t)\}$.
3. Вычислить $v\; \to \; (t\; +\; \Delta \; t)\; =\; v\; \to \; (t)\; +\; 1\; 2\; (a\; \to \; (t)\; +\; a\; \to \; (t\; +\; \Delta \; t))\; \Delta \; t\; \{\backslash \; displaystyle\; \{\backslash \; vec\; \{v\}\}\; (\; t\; +\; \backslash \; Delta\; t)\; =\; \{\backslash \; vec\; \{v\}\}\; (t)\; +\; \{\backslash \; tfrac\; \{1\}\; \{2\}\}\; \backslash ,\; \{\backslash \; big\; (\}\; \{\backslash \; vec\; \{a\}\}\; (t)\; +\; \{\backslash \; vec\; \{a\; \}\}\; (t\; +\; \backslash \; Delta\; t)\; \{\backslash \; big)\}\; \backslash \; Delta\; t\}$.

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

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

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

. Метод Верле по скорости является частным случаем бета-версии Ньюмарка. метод с $\beta \; =\; 0\; \{\backslash \; displaystyle\; \backslash \; beta\; =\; 0\}$и $\gamma \; =\; 1/2\; \{\backslash \; displaystyle\; \backslash \; gamma\; =\; 1/2\}$.

Алгоритмическое представление

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

1 struct Body 2 {3 Vec3d pos {0.0, 0.0, 0.0}; 4 Vec3d vel {2.0, 0.0, 0.0}; // 2 м / с по оси x 5 Vec3d acc {0.0, 0.0, 0.0}; // первые 6 раз ускорения нет double mass = 1.0; // 1 кг 7 двойное перетаскивание = 0,1; // rho * C * Area - упрощенное перетаскивание для этого примера 8 9 / ** 10 * Обновить pos и vel с помощью интеграции "Velocity Verlet" 11 * @param dt DeltaTime / time step [например: 0,01] 12 * / 13 void update (двойной dt) 14 {15 Vec3d new_pos = pos + vel * dt + acc * (dt * dt * 0.5); 16 Vec3d new_acc = apply_forces (); // требуется, только если ускорение непостоянно 17 Vec3d new_vel = vel + (acc + new_acc) * (dt * 0.5); 18 pos = new_pos; 19 vel = new_vel; 20 acc = new_acc; 21} 22 23 Vec3d apply_forces () const 24 {25 Vec3d grav_acc = Vec3d {0.0, 0.0, -9.81}; // 9,81 м / с ^ 2 вниз по оси Z 26 Vec3d drag_force = 0,5 * drag * (vel * abs (vel)); // D = 0.5 * (rho * C * Area * vel ^ 2) 27 Vec3d drag_acc = drag_force / mass; // a = F / m 28 return grav_acc - drag_acc; 29} 30};

Локальная ошибка положения интегратора Верле составляет $O\; (\Delta \; t\; 4)\; \{\backslash \; displaystyle\; O\; (\backslash \; Delta\; t\; ^\; \{4\})\}$как описано выше, а локальная ошибка скорости равна $O\; (\Delta \; t\; 2)\; \{\backslash \; displaystyle\; O\; (\backslash \; Delta\; t\; ^\; \{2\})\}$.

Общая ошибка положения, напротив, составляет $O\; (\Delta \; t\; 2)\; \{\backslash \; displaystyle\; O\; (\backslash \; Delta\; t\; ^\; \{2\})\}$, а глобальная ошибка скорости равна $O\; (\Delta \; t\; 2)\; \{\backslash \; displaystyle\; O\; (\backslash \; Delta\; t\; ^\; \{2\})\}$. Их можно получить, отметив следующее:

$error\; (x\; (t\; 0\; +\; \Delta \; t))\; =\; O\; (\Delta \; t\; 4)\; \{\backslash \; displaystyle\; \backslash \; mathrm\; \{error\}\; \{\backslash \; bigl\; (\}\; x\; (t\_\; \{0\}\; +\; \backslash \; Delta\; t)\; \{\backslash \; bigr)\}\; =\; O\; (\backslash \; Delta\; t\; ^\; \{4\})\}$

и

$x\; (t\; 0\; +\; 2\; \Delta \; t)\; =\; 2\; x\; (t\; 0\; +\; \Delta \; t)\; -\; x\; (\; t\; 0)\; +\; \Delta \; t\; 2\; x\; ″\; (t\; 0\; +\; \Delta \; t)\; +\; O\; (\Delta \; t\; 4).\; \{\backslash \; displaystyle\; x\; (t\_\; \{0\}\; +2\; \backslash \; Delta\; t)\; =\; 2x\; (t\_\; \{0\}\; +\; \backslash \; Delta\; t)\; -x\; (t\_\; \{0\})\; +\; \backslash \; Delta\; t\; ^\; \{2\}\; x\; \text{'}\text{'}\; (t\_\; \{0\; \}\; +\; \backslash \; Delta\; t)\; +\; O\; (\backslash \; Delta\; t\; ^\; \{4\}).\}$

Следовательно,

$error\; (x\; (t\; 0\; +\; 2\; \Delta \; t))\; =\; 2\; error\; (x\; (t\; 0\; +\; \Delta \; t))\; +\; O\; (\Delta \; t\; 4)\; =\; 3\; O\; (\Delta \; t\; 4).\; \{\backslash \; displaystyle\; \backslash \; mathrm\; \{error\}\; \{\backslash \; bigl\; (\}\; x\; (t\_\; \{0\}\; +2\; \backslash \; Delta\; t)\; \{\backslash \; bigr)\}\; =\; 2\; \backslash \; mathrm\; \{error\}\; \{\backslash \; bigl\; (\}\; x\; (t\_\; \{0\}\; +\; \backslash \; Delta\; t)\; \{\backslash \; bigr)\}\; +\; O\; (\backslash \; Delta\; t\; ^\; \{4\})\; =\; 3\; \backslash ,\; O\; (\backslash \; Delta\; t\; ^\; \{4\}).\}$

Аналогично:

$error\; (x\; (t\; 0\; +\; 3\; \Delta \; t))\; знак\; равно\; 6\; О\; (\Delta \; t\; 4),\; \{\backslash \; displaystyle\; \backslash \; mathrm\; \{error\}\; \{\backslash \; bigl\; (\}\; x\; (t\_\; \{0\}\; +3\; \backslash \; Delta\; t)\; \{\backslash \; bigl)\}\; =\; 6\; \backslash ,\; O\; (\; \backslash \; Delta\; t\; ^\; \{4\}),\}$

$ошибка\; (x\; (t\; 0\; +\; 4\; \Delta \; t))\; =\; 10\; O\; (\Delta \; t\; 4),\; \{\backslash \; displaystyle\; \backslash \; mathrm\; \{error\}\; \{\backslash \; bigl\; (\}\; x\; (t\_\; \{0\}\; +4\; \backslash \; Delta\; t)\; \{\backslash \; bigl)\}\; =\; 10\; \backslash ,\; O\; (\backslash \; Delta\; t\; ^\; \{4\}),\}$

$ошибка\; (x\; (t\; 0\; +\; 5\; \Delta \; t))\; =\; 15\; O\; (\Delta \; t\; 4),\; \{\backslash \; displaystyle\; \backslash \; mathrm\; \{error\}\; \{\backslash \; bigl\; (\}\; x\; (t\_\; \{0\}\; +5\; \backslash \; Delta\; t)\; \{\backslash \; bigl)\}\; =\; 15\; \backslash ,\; O\; (\backslash \; Delta\; t\; ^\; \{4\}),\}$

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

$error\; (x\; (t\; 0\; +\; n\; \Delta \; t))\; =\; n\; (n\; +\; 1)\; 2\; O\; (\; \Delta \; t\; 4).\; \{\backslash \; displaystyle\; \backslash \; mathrm\; \{error\}\; \{\backslash \; bigl\; (\}\; x\; (t\_\; \{0\}\; +\; n\; \backslash \; Delta\; t)\; \{\backslash \; bigr)\}\; =\; \{\backslash \; frac\; \{n\; (n\; +\; 1)\}\; \{2\}\}\; \backslash ,\; O\; (\; \backslash \; Delta\; t\; ^\; \{4\}).\}$

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

$ошибка\; (x\; (t\; 0\; +\; T))\; =\; (T\; 2\; 2\; \Delta \; t\; 2\; +\; T\; 2\; \Delta \; t)\; O\; (\Delta \; t\; 4),\; \{\backslash \; displaystyle\; \backslash \; mathrm\; \{error\}\; \{\backslash \; bigl\; (\}\; x\; (t\_\; \{0\}\; +\; T)\; \{\backslash \; bigr)\}\; =\; \backslash \; left\; (\{\backslash \; frac\; \{T\; ^\; \{2\}\}\; \{2\; \backslash \; Delta\; t\; ^\; \{2\}\}\}\; +\; \{\backslash \; frac\; \{T\}\; \{2\; \backslash \; Delta\; t\}\}\; \backslash \; right)\; O\; (\backslash \; Delta\; t\; ^\; \{4\}),\}$

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

$error\; (x\; (t\; 0\; +\; T))\; =\; O\; (\Delta \; t\; 2).\; \{\backslash \; displaystyle\; \backslash \; mathrm\; \{error\}\; \{\backslash \; bigr\; (\}\; x\; (t\_\; \{0\}\; +\; T)\; \{\backslash \; bigl)\}\; =\; O\; (\backslash \; Delta\; t\; ^\; \{2\}).\}$

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

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

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

В одном измерении соотношение между неограниченными позициями $x\; ~\; i\; (t)\; \{\backslash \; displaystyle\; \{\backslash \; tilde\; \{x\}\}\; \_\; \{i\}\; ^\; \{(t)\}\}$и фактические позиции $xi\; (t)\; \{\backslash \; displaystyle\; x\_\; \{i\}\; ^\; \{(t)\}\}$точек i в момент времени t могут быть найдены с помощью алгоритма

$d\; 1\; знак\; равно\; x\; 2\; (t)\; -\; x\; 1\; (t),\; \{\backslash \; displaystyle\; d\_\; \{1\}\; =\; x\_\; \{2\}\; ^\; \{(t)\}\; -\; x\_\; \{1\}\; ^\; \{(t)\},\}$

$d\; 2\; =\; \Vert \; D\; 1\; \Vert ,\; \{\backslash \; displaystyle\; d\_\; \{2\}\; =\; \backslash \; |\; d\_\; \{1\}\; \backslash \; |,\}$

$d\; 3\; =\; d\; 2\; -\; rd\; 2,\; \{\backslash \; displaystyle\; d\_\; \{3\}\; =\; \{\backslash \; frac\; \{d\_\; \{2)\; \}\; -r\}\; \{d\_\; \{2\}\}\},\}$

$x\; 1\; (t\; +\; \Delta \; t)\; =\; x\; ~\; 1\; (t\; +\; \Delta \; t)\; +\; 1\; 2\; d\; 1\; d\; 3,\; \{\backslash \; displaystyle\; x\_\; \{1\}\; ^\; \{(t\; +\; \backslash \; Delta\; t)\}\; =\; \{\backslash \; tilde\; \{x\}\}\; \_\; \{1\}\; ^\; \{(t\; +\; \backslash \; Delta\; t)\}\; +\; \{\backslash \; frac\; \{1\}\; \{2\}\}\; d\_\; \{1\}\; d\_\; \{3\},\}$

$х\; 2\; (t\; +\; \Delta \; t)\; =\; x\; ~\; 2\; (t\; +\; \Delta \; t)\; -\; 1\; 2\; d\; 1\; d\; 3.\; \{\backslash \; displaystyle\; x\_\; \{2\}\; ^\; \{(t\; +\; \backslash \; Delta\; t)\}\; =\; \{\backslash \; tilde\; \{x\}\}\; \_\; \{2\}\; ^\; \{(t\; +\; \backslash \; Delta\; t)\}\; -\; \{\backslash \; frac\; \{1\}\; \{2\}\}\; d\_\; \{\; 1\}\; d\_\; \{3\}.\}$

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

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

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

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

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

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

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

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

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

1. ^Верле, Лу (1967). «Компьютерные« эксперименты »над классическими жидкостями. I. Термодинамические свойства молекул Леннарда-Джонса». Физический обзор. 159 (1): 98–103. Bibcode : 1967PhRv..159... 98V. doi : 10.1103 / PhysRev.159.98.
2. ^Press, W.H.; Теукольский, С. А.; Vetterling, W. T.; Фланнери, Б. П. (2007). «Раздел 17.4. Консервативные уравнения второго порядка». Числовые рецепты: искусство научных вычислений (3-е изд.). Нью-Йорк: Издательство Кембриджского университета. ISBN 978-0-521-88068-8.
3. ^веб-страница Архивировано 2004-08-03 на Wayback Machine с описание метода Штёрмера.
4. ^Даммер, Джонатан. «Простой метод интеграции Верле с временной коррекцией».
5. ^Свуп, Уильям К.; Х. К. Андерсен; П. Х. Беренс; К. Р. Уилсон (1 января 1982 г.). «Метод компьютерного моделирования для расчета констант равновесия для образования физических кластеров молекул: приложение к небольшим кластерам воды». Журнал химической физики. 76 (1): 648 (Приложение). Bibcode : 1982JChPh..76..637S. doi : 10.1063 / 1.442716.
6. ^Хайрер, Эрнст; Любич, Кристиан; Ваннер, Герхард (2003). «Геометрическое численное интегрирование, иллюстрированное методом Штёрмера / Верле». Acta Numerica. 12 : 399–450. Bibcode : 2003AcNum..12..399H. CiteSeerX 10.1.1.7.7106. doi : 10.1017 / S0962492902000144.
7. ^Руководство пользователя SuperLU.
8. ^Baraff, D.; Виткин, А. (1998). «Большие шаги в моделировании ткани» (PDF). Компьютерная графика Труды. Ежегодная серия конференций: 43–54.

• Демонстрация интеграции Verlet и код в виде Java-апплета
• Расширенная физика персонажей Томасом Якобсеном
• Теория моделирования молекулярной динамики - внизу страницы