Совершенствование алгоритмов численного моделирования в методе динамики частиц

Задача

• ...

Построение модифицированного метода Рунге-Кутты 4 порядка

Рассмотрим задачу Коши

[math] dx/dt=f(x,t) \ \ (1) [/math]

[math] x(0)=x_0 [/math]

Для неизвестной вектор-функции x(t), в качестве которой для примера может быть взят вектор [math] (r,v) = (x,y,z,u,v,w) [/math] координат позиции и скорости тела. Данная задача может быть решена численно классическим методом Рунге-Кутты четвёртого порядка.

[math] k_1=\Delta t f(x_n,t_n ) [/math]

[math] k_2=\Delta t f(x_n+k_1/2,t_n+Δt/2) [/math]

[math] k_3=\Delta t f(x_n+k_2/2,t_n+Δt/2) [/math]

[math] k_4=\Delta t f(x_n+k_3,t_n+Δt) [/math]

[math] x_{n+1}=x_n+1/6 (k_1+2k_2+2k_3+k_4 ) \ \ (2) [/math]

По сравнению с методами Эйлера, Лагранжа и Верле, данный метод имеет более высокий порядок точности. Однако классический метод Рунге-Кутты четвёртого порядка имеет одну особенность, связанную с необходимостью вычислять функцию [math] f(x,t) [/math] четыре раза за одну временную итерацию. Потому этот метод становится неэффективным в вычислительных задачах, где основное расчётное время тратится на вычисление правой части системы дифференциальных уравнений, как, например, это имеет место в случае расчёта молекулярно-динамической задачи множества частиц. Вследствие данной особенности применение метода Рунге-Кутты становится неэффективным и даже его исключительная точность теряет свою значимость.

Ниже приводится модификация метода Рунге-Кутты 4 порядка, где с помощью одного хитрого приёма удаётся избежать многократного вычисления функцию [math] f(x,t) [/math] на одном временном шаге и в то же время сохранить высокий порядок по времени.

Идея заключается в разложении функций [math] f(x_n+k_i/2,t_n+Δt/2) [/math] в ряд Тейлора в окрестности точки [math] (x_n,t_n) [/math].

[math] f(x_n+k_i/2,t_n+Δt/2)= \frac {\partial f} {\partial x} (x_n,t_n )∙k_i/2+ \frac {\partial f} {\partial t} (x_n,t_n ) Δt/2 + ... \ \ (3) [/math]

Здесь присутствуют малоприятные производные, однако, как потом окажется, с ними можно будет легко разобраться. Сколько членов в разложении нужно оставить, чтобы в схеме сохранился четвёртый порядок? – До [math] (Δt)^4 [/math] и [math] (k_i )^4 [/math] или меньше?

Для слагаемых с локальными производными по времени ответ очевиден – необходимо удерживать всё вплоть до [math] (Δt)^4 [/math], ибо в противном случае мы потеряем наш 4-й порядок по времени для схемы в целом. Однако для [math] k_i [/math] на самом деле достаточно только первой производной.

В случае, когда правая часть (1) не зависит явно от времени, (3) предельно упрощается.

[math] f(x_n+k_i/2,t_n+Δt/2)=∂f/∂x (x_n,t_n )∙k_i/2 \ \ (4) [/math]

Данная ситуация имеет место при молекулярно-динамическом моделировании, поскольку потенциал взаимодействия, как правило, является функцией только координат и скоростей частиц.

Запишем (4) в случае молекулярно-динамического моделирования. В нашем случае неизвестная вектор-функция [math] x = {r,v} [/math].

[math] (x,t)=(v,\frac{F(r)}{m}) [/math]

[math] f(x_n+k_i/2,t_n+Δt/2) = (v,\frac{F(r+k_i^r/2)}{m}) = ( v+\frac {k_1^v}{2},\frac {F(r+\frac {k_1^r}{2})}{m} )= [/math]

[math] =( v+\frac {k_1^v}{2},\frac {F_n + \frac {dF}{dr} (r) \frac{k_1^r}{2}}{m} ) \ \ (5) [/math]

где

[math] \frac {dF}{dr} = \frac {\Delta F}{\Delta r} = \frac {\Delta F}{\Delta r} \frac {\Delta t}{\Delta t} = \frac {\Delta F}{\Delta t} \frac {\Delta t}{\Delta r} =\frac {F_n - F_{n-1}}{\Delta t } \frac {1}{v} \ \ (6) [/math]

где [math]F_n[/math] и [math]F_{n-1}[/math] - силы, действующие на рассматриваемую частицу на данном временном слое [math]n[/math] и на предыдущем временном слое [math](n-1)[/math] соответственно. Это значит, что в памяти необходимо хранить силы с прошлого временного слоя, что не очень хорошо. Это главный, но не очень существенный минус данного метода. С учётом вышесказанного, мы можем записать конечно-разностную схему Рунге-Кутта для молекулярно-динамического моделирования

[math]k_1^v=\Delta t \frac {F_n}{m}[/math]

[math]k_2^v=\Delta t \frac {F_n}{m}+\frac{F_n-F_(n-1)}{2m}[/math]

[math]k_3^v=\Delta t \frac {F_n}{m}+\frac{F_n-F_(n-1)}{2m}[/math]

[math]k_4^v=\Delta t \frac {F_n}{m}+\frac{F_n-F_(n-1)}{m}[/math]

[math]k_1^r=\Delta t \frac {F_n}{m}[/math]

[math]k_2^r=\Delta tv_n+\frac {v_n-v_(n-1)}{2}=\Delta t \frac {v_n+(k_1^v}{2}[/math]

[math]k_3^r=\Delta tv_n+\frac {v_n-v_(n-1)}{2}= \Delta t \frac {v_n+(k_2^v}{2})[/math]

[math]k_4^r=\Delta tv_n+v_n-v_(n-1)= \Delta t(v_n+k_3^v )[/math]

[math]r_{n+1}=r_n+1/6 (k_1^r+2k_2^r+2k_3^r+k_4^r )[/math]

[math]v_{n+1}=v_n+1/6 (k_1^((v) )+2k_2^((v) )+2k_3^((v) )+k_4^((v) ) ) \ \ (8) [/math]

Выражения для [math] k_i^r[/math] представлены в двух видах: один – с учётом рассмотренного упрощения, а второй - классический Рунге-Куттовский, который применим для уравнения по [math] r [/math] ввиду простой правой части.

Обезразмеривание системы как способ уменьшения накопления вычислительной ошибки

• ...

Frozen Particles & Press Particles

• ...