Редактирование: Совершенствование алгоритмов численного моделирования в методе динамики частиц
Внимание! Вы не авторизовались на сайте. Ваш IP-адрес будет публично видимым, если вы будете вносить любые правки. Если вы войдёте или создадите учётную запись, правки вместо этого будут связаны с вашим именем пользователя, а также у вас появятся другие преимущества.
Правка может быть отменена. Пожалуйста, просмотрите сравнение версий, чтобы убедиться, что это именно те изменения, которые вас интересуют, и нажмите «Записать страницу», чтобы изменения вступили в силу.
Текущая версия | Ваш текст | ||
Строка 1: | Строка 1: | ||
− | == | + | == Задача == |
− | * | + | * ... |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
== Построение модифицированного метода Рунге-Кутты 4 порядка == | == Построение модифицированного метода Рунге-Кутты 4 порядка == | ||
Строка 27: | Строка 15: | ||
<math> k_1=\Delta t f(x_n,t_n ) </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+ | + | <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+ | + | <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+ | + | <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> x_{n+1}=x_n+1/6 (k_1+2k_2+2k_3+k_4 ) \ \ (2) </math> | ||
Строка 39: | Строка 27: | ||
Ниже приводится модификация метода Рунге-Кутты 4 порядка, где с помощью одного хитрого приёма удаётся избежать многократного вычисления функцию <math> f(x,t) </math> на одном временном шаге и в то же время сохранить высокий порядок по времени. | Ниже приводится модификация метода Рунге-Кутты 4 порядка, где с помощью одного хитрого приёма удаётся избежать многократного вычисления функцию <math> f(x,t) </math> на одном временном шаге и в то же время сохранить высокий порядок по времени. | ||
− | Идея заключается в разложении функций <math> f(x_n+k_i/2,t_n+ | + | Идея заключается в разложении функций <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+ | + | <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> ( | + | Здесь присутствуют малоприятные производные, однако, как потом окажется, с ними можно будет легко разобраться. Сколько членов в разложении нужно оставить, чтобы в схеме сохранился четвёртый порядок? – До <math> (Δt)^4 </math> и <math> (k_i )^4 </math> или меньше? |
− | Для слагаемых с локальными производными по времени ответ очевиден – необходимо удерживать всё вплоть до <math> ( | + | Для слагаемых с локальными производными по времени ответ очевиден – необходимо удерживать всё вплоть до <math> (Δt)^4 </math>, ибо в противном случае мы потеряем наш 4-й порядок по времени для схемы в целом. Однако для <math> k_i </math> на самом деле достаточно только первой производной. |
В случае, когда правая часть (1) не зависит явно от времени, (3) предельно упрощается. | В случае, когда правая часть (1) не зависит явно от времени, (3) предельно упрощается. | ||
− | <math> f(x_n+k_i/2,t_n+ | + | <math> f(x_n+k_i/2,t_n+Δt/2)=∂f/∂x (x_n,t_n )∙k_i/2 \ \ (4) </math> |
Данная ситуация имеет место при молекулярно-динамическом моделировании, поскольку потенциал взаимодействия, как правило, является функцией только координат и скоростей частиц. | Данная ситуация имеет место при молекулярно-динамическом моделировании, поскольку потенциал взаимодействия, как правило, является функцией только координат и скоростей частиц. | ||
Строка 57: | Строка 45: | ||
<math> (x,t)=(v,\frac{F(r)}{m}) </math> | <math> (x,t)=(v,\frac{F(r)}{m}) </math> | ||
− | <math> f(x_n+k_i/2,t_n+ | + | <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> =( v+\frac {k_1^v}{2},\frac {F_n + \frac {dF}{dr} (r) \frac{k_1^r}{2}}{m} ) \ \ (5) </math> | ||
Строка 93: | Строка 81: | ||
Существует четыре источника погрешностей результата численного метода: | Существует четыре источника погрешностей результата численного метода: | ||
− | + | ** погрешность математической модели, | |
− | + | ** погрешность исходных данных, | |
− | + | ** погрешность метода (в литературе ее также называют погрешностью обрыва), | |
− | + | ** погрешность округления (машинная погрешность). | |
− | |||
− | |||
− | |||
− | |||
− | |||
Каждый из этих источников по-своему специфичен и интересен, однако нас здесь будет волновать последний, четвёртый источник, связанный с несовершенством вычислительного оборудования. Главные причины больших случайных погрешностей ЭВМ: | Каждый из этих источников по-своему специфичен и интересен, однако нас здесь будет волновать последний, четвёртый источник, связанный с несовершенством вычислительного оборудования. Главные причины больших случайных погрешностей ЭВМ: | ||
− | + | ** метод обрыва и округления, принятый в машине; | |
− | + | ** потеря значащих разрядов при вычитании; | |
− | + | ** техническое состояние машины; | |
− | + | ** потеря разрядов при превышении допустимой разрядности представления чисел (например, при делении на маленькие числа). | |
− | |||
− | |||
− | |||
− | |||
− | |||
Проблема ошибок такого рода решается как правило заменой ЭВМ на более мощную. Однако это не единственный способ решения проблемы. Существует ряд алгоритмических приёмов, позволяющих добиться меньшей машинной погрешности, работая на прежней ЭВМ. Один из наиболее эффективных приёмов - обезразмеривание. | Проблема ошибок такого рода решается как правило заменой ЭВМ на более мощную. Однако это не единственный способ решения проблемы. Существует ряд алгоритмических приёмов, позволяющих добиться меньшей машинной погрешности, работая на прежней ЭВМ. Один из наиболее эффективных приёмов - обезразмеривание. | ||
− | Обезразмеривание как таковое возникло намного раньше ЭВМ и вовсе не из желания сократить машинную ошибку, о которой тогда и знать не знали. Первоначально обезразмеривание было исключительно математическим приёмом, позволяющим перейти от физической задачи (в которой присутствовали физические величины, обладающие смыслом и размерностью) к абстрактной математической задаче (в которой величины - просто некоторые отвлечённые параметры функции), что | + | Обезразмеривание как таковое возникло намного раньше ЭВМ и вовсе не из желания сократить машинную ошибку, о которой тогда и знать не знали. Первоначально обезразмеривание было исключительно математическим приёмом, позволяющим перейти от физической задачи (в которой присутствовали физические величины, обладающие смыслом и размерностью) к абстрактной математической задаче (в которой величины - просто некоторые отвлечённые параметры функции), что давао возможность применять формальные математическое методы решения уравнений, например метод малого параметра. |
Обезразмеривание позволяет перейти от физических величин, которые могут иметь очень большие (модуль Юнга) или очень малые (постоянная Планка) по модулю значения к их безразмерным аналогам, которые порядка единиц и в крайнем случае десятков. Понятно что второй случай более предпочтителен в вычислительном плане, ибо при оперировании данными безразмерными величинами происходит гораздо более медленное накопление ошибки. | Обезразмеривание позволяет перейти от физических величин, которые могут иметь очень большие (модуль Юнга) или очень малые (постоянная Планка) по модулю значения к их безразмерным аналогам, которые порядка единиц и в крайнем случае десятков. Понятно что второй случай более предпочтителен в вычислительном плане, ибо при оперировании данными безразмерными величинами происходит гораздо более медленное накопление ошибки. | ||
Строка 122: | Строка 100: | ||
Рассмотрим задачу Коши | Рассмотрим задачу Коши | ||
− | <math> \frac{dv_i}{dt}=\frac{f_i ( | + | <math> \frac{dv_i}{dt}=\frac{f_i (r_ij,v_ij )}{m_i} </math> |
<math> \frac{dr_i}{dt}=v_i </math> | <math> \frac{dr_i}{dt}=v_i </math> | ||
Строка 136: | Строка 114: | ||
Сделаем замену переменных | Сделаем замену переменных | ||
− | <math>r_i | + | <math>r_i -> LR_i</math> |
− | <math> v_i | + | <math> v_i -> L/τ V_i</math> |
− | <math>t | + | <math>t -> τ </math> |
− | <math> | + | <math>f_ -> \m frac{L}{τ^2} F_i </math> |
<math>i=1,..,N</math> | <math>i=1,..,N</math> | ||
где <math>\tau, R_i,V_i,F_i </math> – новые безразмерные параметры (время, радиус-вектор, вектор скорости и сила соответственно). | где <math>\tau, R_i,V_i,F_i </math> – новые безразмерные параметры (время, радиус-вектор, вектор скорости и сила соответственно). | ||
− | + | Система (1) в безразмерном виде будет таковой | |
− | + | (dV_i)/dτ=F_i (R_ij,V_ij ) | |
− | + | (dR_i)/dτ=V_i | |
− | + | R_i (0)=(r_i^0)/L | |
− | + | V_i (0)=(v_i^0 τ)/L | |
− | + | i=1,..,N (3) | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
Однако на этом обезразмеривание ещё не заканчивается. | Однако на этом обезразмеривание ещё не заканчивается. | ||
− | + | f_i (r_ij )=∑_(j=1)^N▒〖f(r_ij)〗 (4) | |
− | |||
− | |||
где | где | ||
− | + | f(r_ij )=-∂U/(∂r_ij ) (5) | |
− | + | Если в качестве потенциала взаимодействия U взять потенциал Леннарда-Джонса | |
− | + | U(r_ij )=D[(a/r_ij )^12-2(a/r_ij )^6 ]=D[(1/R_ij )^12-2(1/R_ij )^6 ] (6) | |
− | Если в качестве потенциала взаимодействия U взять потенциал | + | Глубина D потенциальной ямы имеет размерность энергии. Если мы хотим обезразмерить (6), то нам достаточно ввести новую безразмерную глубину D ̃ потенциальной ямы по правилу |
− | + | D=D ̃ (mL^2)/τ^2 (7) | |
− | |||
− | |||
− | Глубина | ||
− | |||
− | |||
− | |||
Тогда новая обезразмеренная потенциальная энергия запишется в виде | Тогда новая обезразмеренная потенциальная энергия запишется в виде | ||
− | + | U ̃(R_ij )=D ̃[(1/R_ij )^12-2(1/R_ij )^6 ] (8) | |
− | + | А безразмерная сила F ̃ будет выражаться через потенциал привычной формулой | |
− | + | F(R_ij )=-(∂U ̃)/(∂R_ij ) (9) | |
− | А безразмерная сила | ||
− | |||
− | |||
== Frozen Particles & Press Particles == | == Frozen Particles & Press Particles == | ||
− | + | * ... | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− |