Численный алгоритм решения системы уравнений движения

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

Рассмотрим данный алгоритм на примере одномерного движения частицы. В двумерном и трехмерном случае появляются еще уравнения для y, v_y,a_y и z, v_z, a_z.

пусть x_n=x(t_n), x_n+1=x(t_n+1) и т.д.

По формуле Тейлора имеем для x_n+1 и v_n+1

x_n+1=x(t_n+Δt)= x_n + x́(t_n)Δt+1/2 x́́́ ́(t_n)(Δt)²+O((Δt)³) (1);

v_n+1= v_n+v́(t_n) Δt+ O((Δt)²) (2).

Учтем, что x́(t_n)= v_n и x́́́ ́(t_n)= v́(t_n)=a_n, тогда

x_n+1=x(t_n+Δt)= x_n + v_nΔt+1/2 a_n(Δt)²+O((Δt)³) (3);

v_n+1= v_n+ a_nΔt+ O((Δt)²) (4).

Если использовать формулы (3) и (4) будет большая ошибка. Рассмотрим

x_n-1=x(t_n-Δt)≈ x_n - v_nΔt+1/2 a_n(Δt)² (5).

Сложим левые и правые части (3) и (5) и выразим x_n+1 из получившегося равенства, получим

x_n+1+ x_n-1≈2 x_n+ a_n(Δt)²,

x_n+1≈2 x_n- x_n-1+ a_n(Δt)² (6).

Вычтем из (3) равенство (5) и выразим из получившегося равенства v_n:

v_n= (7).

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

Рассмотрим равенство (6). В правой части добавим и вычтем 0,5x_n+1:

x_n+1≈2 x_n- x_n-1+ a_n(Δt)²+0,5x_n+1-0,5x_n+1= x_n+0,5(x_n+1- x_n-1)-0,5 x_n-1-0,5 x_n+1+ x_n+ a_n(Δt)²=

=x_n+ v_nΔt-0,5(x_n+1-2 x_n + x_n-1) + a_n(Δt)² (8).

Из (6) выразим a_n

a_n= (9).

Поэтому (8) с учетом (9) имеет вид

x_n+1= x_n+ v_nΔt+0,5a_n(Δt)² (10).

Аналогично из (7)

(11)

из (10) =x_n+1+ v_n+1Δt+0,5a_n+1(Δt)² подставим в (11). Получим

Выразив из последнего равенства v_n+1, получим

v_n+1= v_n+0,5(a_n+1 +a_n)Δt (12)

(10) и (12) – скоростная форма алгоритма Верле. Она является самостартующей.

Таким образом, алгоритм выглядит следующим образом:

· Находят x_n+1= x_n+ v_nΔt+0,5a_n(Δt)²;

· Вычисляют ускорение a_n+1 из СДУ, используя в правой части уравнений данные с x_n+1;

· По найденному a_n+1 подправляют скорость

v_n+1= v_n+0,5(a_n+1 +a_n)Δt

Дата добавления: 2014-01-06