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

Задача

• ...

Построение модифицированного метода Рунге-Кутты 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 ( v_n+\frac{k_1^v}{2}) [/math]

[math]k_3^r=\Delta tv_n+\frac {v_n-v_{n-1}}{2}= \Delta t ( v_n+\frac{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] ввиду простой правой части.

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

Существует четыре источника погрешностей результата численного метода: 1) погрешность математической модели, 2) погрешность исходных данных, 3) погрешность метода (в литературе ее также называют погрешностью обрыва), 4) погрешность округления (машинная погрешность).

Каждый из этих источников по-своему специфичен и интересен, однако нас здесь будет волновать последний, четвёртый источник, связанный с несовершенством вычислительного оборудования. Главные причины больших случайных погрешностей ЭВМ: 1) метод обрыва и округления, принятый в машине; 2) потеря значащих разрядов при вычитании; 3) техническое состояние машины; 4) потеря разрядов при превышении допустимой разрядности представления чисел (например, при делении на маленькие числа).

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

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

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

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

[math] \frac{dv_i}{dt}=\frac{f_i (r_{ij},v_{ij} )}{m_i} [/math]

[math] \frac{dr_i}{dt}=v_i [/math]

[math] r_i (0)=r_i^0 [/math]

[math]v_i (0)=v_i^0 [/math]

[math]i=1,..,N[/math]

Пусть [math]\tau[/math] – характерное время протекания рассматриваемых процессов, [math]a[/math] – характерная длина расчётной области, [math]m[/math] – характерная масса частиц.

Сделаем замену переменных

[math]r_i = a R_i[/math]

[math] v_i = \frac{a}{\tau} V_i[/math]

[math]t = \tau [/math]

[math] f_i = m \frac{a}{τ^2} F_i [/math]

[math]i=1,..,N[/math]

где [math]\tau, R_i,V_i,F_i [/math] – новые безразмерные параметры (время, радиус-вектор, вектор скорости и сила соответственно).

Исходная система Коши в безразмерном виде будет таковой

[math]\frac{dV_i}{d\tau}=F_i (R_{ij},V_{ij} )[/math]

[math]\frac{dR_i}{d\tau}=V_i[/math]

[math] R_i (0)=(r_i^0)/L[/math]

[math]V_i (0)=(v_i^0 \tau)/L[/math]

[math]i=1,..,N[/math]

Однако на этом обезразмеривание ещё не заканчивается.

[math] f_i (r_{ij} )=\sum_{j=0}^{N}f(r_{ij})[/math]

где

[math] f(r_{ij} )=- \frac{\partial u}{\partial r_{ij}} [/math]

Если в качестве потенциала взаимодействия U взять потенциал Леннарда-Джонса

[math] u(r_{ij} )=d[(\frac{a}{r_{ij}} )^{12}-2(\frac{a}{r_{ij}} )^6 ]=D[(\frac{1}{R_{ij}} )^{12}-2( \frac{1}{R_{ij}} )^6 ] [/math]

Глубина [math]d[/math] потенциальной ямы имеет размерность энергии. Если мы хотим обезразмерить потенциал, то нам достаточно ввести новую безразмерную глубину [math] D [/math] потенциальной ямы по правилу

[math]d=D \frac{mL^2}{τ^2}[/math]

Тогда новая обезразмеренная потенциальная энергия запишется в виде

[math]U(R_ij )=D ((\frac{1}{R_{ij}} )^{12}-2(\frac{1}{R_{ij}} )^6 ) [/math]

А безразмерная сила [math]F[/math] будет выражаться через потенциал привычной формулой

[math]F(R_{ij} )- \frac{\partial U}{\partial R_{ij}} [/math]

Frozen Particles & Press Particles

• ...