Курс: Общий физический практикум Сегодня: _________________ 2009 г. Склярова Елена Александровна

Лекция 5 Тема: Численное моделирование 1. Численное моделирование 2. Численное интегрирование уравнений Ньютона (алгоритм Верле, Бимана, Шофилда, Рунге-Кутты) 3. Динамика систем многих частиц Содержание лекции: Сегодня: __________________ 2009 г.

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

Численное интегрирование уравнений Ньютона Рассмотрим одномерное движение частицы и запишем уравнения Ньютона в виде (1) (2) Целью всех конечно-разностных методов является вычисление значений x n+1 и n+1 (точка в «фазовом пространстве») в момент времени t n+1 = t n + t.

Численное интегрирование уравнений Ньютона !!! Величину шага t надо выбирать таким образом, чтобы метод интегрирования порождал устойчивое решение. 1 способ проверки устойчивости метода заключается в контроле величины полной энергии и обеспечении того, чтобы она не отклонялась от начального значения в случае, когда полная энергия сохраняется. Достаточно большое значение шага t приводит к несохранению полной энергии и неустойчивым решениям для x n+1 и n+1, т.е. к таким численным решениям, которые все больше отклоняются с течением времени от истинного решения.

Суть многих алгоритмов можно понять, разлагая v n+1 v(t n + t) и x n+1 x(t n + t) в ряд Тейлора. Запишем v n+1 = v n +a n t + O[( t) 2 ](3) и (4) Хорошо известный метод Эйлера эквивалентен сохранению в (3-4) членов О( t): v n+1 = v n + a n t(5) и (6) Численное интегрирование уравнений Ньютона

Поскольку мы удержали в формулах (5-6) члены порядка t, то «локальная» погрешность (погрешность на шаге) составляет величину O( t) 2. Однако от шага к шагу погрешности накапливаются, поэтому можно предполагать, что «глобальная» погрешность, представляющая собой суммарную погрешность за рассматриваемый промежуток времени, будет величиной O( t). Эта оценка погрешности вполне правдоподобна, поскольку число шагов, на которое разбивается временной интервал, пропорционально 1/ t. Следовательно, порядок глобальной погрешности увеличивается в t раз по отношению к локальной погрешности. Поскольку принято говорить, что метод имеет n-й порядок аппроксимации, если его локальная погрешность равна O(( t)n+1), то метод Эйлера относится к методам первого порядка. Численное интегрирование уравнений Ньютона

Метод Эйлера является асимметричным, поскольку он продвигает решение на один временной шаг t, а использует при этом информацию о производной только в начальной точке интервала! Для полноты изложения повторим алгоритм Эйлера- Кромера или приближение по «последней точке»: v n+1 = v n + a n t(7) и (8) Численное интегрирование уравнений Ньютона

Пожалуй, самый очевидный путь улучшения метода Эйлера состоит в использовании для вычисления нового значения координаты средней на отрезке скорости. Соответствующий метод средней точки можно записать в виде v n+1 = v n + a n t(9) и (10) Заметим, что если подставить выражение (9) для v n+1 в (10), то получим (11)

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

Численное интегрирование уравнений Ньютона Метод полушага относится к методам более высокого порядка точности с ограниченной погрешностью. Принимается, что средняя скорость на отрезке равна значению скорости в середине отрезка. Метод полушага можно записать в виде v n+1/2 = v n-1/2 + a n Δt (12 а) и x n+1 = x n + v n+1/2 Δt.(12 б)

Численное интегрирование уравнений Ньютона Заметим, что метод полушага не является «самостартующим», т.е. формулы (12) не позволяют вычислить v 1/2. Это неудобство можно преодолеть, положив (12 в) Поскольку формулы (12) можно повторять до бесконечности, то метод полушага получил широкое распространение в учебной литературе

Один из наиболее известных алгоритмов высокого порядка приписывается Верле. По аналогии с (3-4) запишем разложение в ряд Тейлора для x n-1 : (13) Если сложить формулы интегрирования вперед и назад выражения (3-4) и (12) соответственно, то получим (14) Или (15 а) Численное интегрирование уравнений Ньютона

Аналогично вычитание разложений в ряд Тейлора для x n+1 и x n-1 дает (15 б) При этом связанная с алгоритмом Верле (15) глобальная погрешность имеет третий порядок для координаты и второй порядок для скорости. Однако скорость не участвует в интегрировании уравнений движения. В литературе по численному анализу алгоритм Верле называется «неявной симметричной разностной схемой».

Численное интегрирование уравнений Ньютона Менее известной, но математически эквивалентной версией алгоритма Верле является схема (16 а) и (16 б) Видно, что схема (16), называемая скоростной формой алгоритма Верле, является самостартующей и не приводит к накоплению погрешности округления.

Численное интегрирование уравнений Ньютона Формулы (16) можно «вывести» из (15) следующим образом. Сначала прибавим и вычтем из обеих частей равенства (15а) по (1/2)x n+1 и запишем (17)

Численное интегрирование уравнений Ньютона Из (15а) найдем a n для метода Верле: (18) Легко видеть, что подстановка (18) в выражение (17) приводит к (16 а). В том же духе перепишем (15а) для v n+1 : (19)

Теперь перепишем формулу (15 а) для x n+2 и подставим полученный результат в (19). Имеем (20) Затем, используя выражение (15а) для x n+1, повторим эту процедуру и подставим x n+1 в (20); после несложных выкладок получаем требуемый результат (16б). Численное интегрирование уравнений Ньютона

Другой полезный алгоритм, в котором нет накопления погрешности округления, присуще алгоритму Верле, принадлежит Биману и Шофилду. Запишем алгоритм Бимана в следующем виде: (21а) и (21б) Заметим, что точность расчета траекторий по схеме (21) не выше, чем в алгоритме Верле. Ее преимущество заключается, пожалуй, в том, что обычно она лучше сохраняет энергию.

Численное интегрирование уравнений Ньютона Завершим наше обсуждение кратким изложением двух методов, которые обычно приводятся в учебниках по численному анализу.

Численное интегрирование уравнений Ньютона Один пример метода предиктор-корректор определяется следующим образом. Сначала мы предсказываем (predictor) новое значение координаты: предиктор: (22 а) Предсказанное значение координаты позволяет определить ускорение. Затем, используя, получаем скорректированные (corrector) значения v n+1 и x n+1 : корректор: (22 б)

Для объяснения метода Рунге-Кутты рассмотрим сначала решение дифференциального уравнения первого порядка (23) Метод Рунге-Кутты второго порядка для решения уравнения (23) можно, пользуясь стандартными обозначениями, записать в следующем виде: (24) Численное интегрирование уравнений Ньютона

Смысл формул (24) состоит в следующем. В методе Эйлера предполагается, что для экстраполяции в следующую точку можно использовать наклон кривой f(x n, y n ) в точке (x n, y n ), т.е. y n+1 = y n + f(x n, y n ) x. Однако можно повысить точность оценки наклона, если методом Эйлера провести экстраполяцию в срединную точку отрезка, а затем использовать центральную производную на всем отрезке. Отсюда оценка наклона в методе Рунге-Кутты равна где

Численное интегрирование уравнений Ньютона Применение метода Рунге-Кутты к уравнениям движения Ньютона (1) дает (25)

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

Этот пример с чашкой кофе ставит нас перед проблемой: как можно, исходя из известных межмолекулярных взаимодействий, понять наблюдаемое поведение сложной многочастичной системы? Самый очевидный подход – решить эту задачу в лоб, моделируя на компьютере саму задачу многих частиц. На самом деле этот подход, называемый методом молекулярной динамики, применили к «небольшим» системам многих частиц, насчитывающим обычно от нескольких сотен до нескольких тысяч частиц, и он уже много дал для понимания наблюдаемых свойств газов, жидкостей и твердых тел. Однако детальное знание траекторий 10 4 или даже частиц ничего не даст, если мы не знаем, какие именно вопросы требуют выяснения. Какие основные свойства и закономерности проявляют системы многих тем?

Потенциал межмолекулярного взаимодействия Первым делом нам необходимо определить модель системы, которую мы желаем моделировать. Поскольку мы хотим понять качественные свойства систем многих частиц, пойдем на упрощение задачи, предполагая, что динамику можно считать классической, а молекулы – химически инертными шариками. Мы предполагаем также, что сила взаимодействия любых двух молекул зависит только от расстояния между ними. В этом случае полная потенциальная энергия U определяется суммой двухчастичных взаимодействий: (1) где V(r ij ) зависит только от абсолютной величины расстояния r ij между частицами i и j. Парное взаимодействие вида (1) соответствует «простым» жидкостям, например жидкому аргону.

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

Одной из наиболее употребительных феноменологических формул для V(r) является потенциал Леннарда-Джонса (2) График потенциала Леннарда-Джонса

Хотя зависимость r -6 в формуле (2) получена теоретически, зависимость r -12 выбрана только из соображений удобства. Потенциал Леннарда-Джонса параметризуется «длиной» и «энергией». Заметим, что при r = V(r) = 0. Параметр представляет собой глубину потенциальной ямы в точке минимума V(r); этот минимум расположен при расстоянии r = 2 1/6. Заметим, что данный потенциал является «короткодействующим» и V(r) для r > 2.5 по существу равен нулю. Потенциал межмолекулярного взаимодействия

Удобно выражать длины, энергию и массу в единицах, и m, где m – масса частицы. Мы измеряем скорости в единицах ( /m) 1/2, а время в – в единицах = (m 2 ) 1/2. Для жидкого аргона параметры и потенциала Леннарда-Джонса составляют /k B = K и = Å. Масса атома аргона равна гр, и отсюда = с. Во избежание путаницы отметим все безразмерные или приведенные величины звездочкой. Например, приведенная двумерная плотность определяется соотношением * = / 2.

Численный алгоритм Теперь, когда у нас есть четко описанная модель системы многих частиц, необходимо познакомиться с каким-нибудь методом численного интегрирования для расчета траекторий каждой частицы. А лгоритм только для одной компоненты движения частицы (алгоритм Верле в скоростной форме)

Поскольку новая координата x n+1 вычисляется с использованием не только скорости v n, но и ускорения a n, алгоритм Верле обладает «более высоким порядком» по t, чем алгоритмы Эйлера и Эйлера-Кромера. Новая координата используется для нахождения нового ускорения a n+1, которое вместе с a n используется для получения новой скорости v n+1. Численный алгоритм

Численный алгоритм. Краевые условия Полезное моделирование должно включать в себя все характерные особенности рассматриваемой физической системы. К онечной целью модельных расчетов является получение оценок поведения макроскопических систем, т.е. систем, содержащих порядка N частиц. Рассмотрим сферический резервуар с водой. Доля молекул воды вблизи стенок пропорциональна отношению поверхности к объему (4 R 2 )/(4 R 3 /3). Поскольку N = (4/3 R 3 ), где - плотность, доля частиц вблизи стенок пропорциональна N 2/3 /N = N -1/3, что при N пренебрежимо мало. По сравнению с этим количество частиц, которое можно изучать в моделях молекулярной динамики, составляет обычно , и доля частиц вблизи стенок не мала.

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

Один из способов минимизировать поверхностные эффекты и более точно промоделировать свойства макроскопической системы заключается в использовании периодических краевых условий. Реализация периодических краевых условий для короткодействующих взаимодействий, таких как потенциал Леннарда-Джонса, хорошо знакома всем играющим в видеоигры. Рассмотрим сначала одномерный «ящик», содержащий N частиц, движение которых ограничено отрезком прямой длиной L. Концы отрезка прямой играют роль «стенок».

Численный алгоритм. Краевые условия Использование периодических краевых условий эквивалентно сворачиванию этого отрезка в кольцо. Заметим, что, поскольку расстоянию между частицами соответствует отрезок дуги, максимальное расстояние равно L/2.

В двумерном случае мы можем представить себе ящик, у которого противоположные ребра соединены так, что ящик превращается в поверхность тора (форма тороидальной камеры). Поэтому, если частица пересекает ребро ящика, она входит заново в противолежащее ребро. При этом максимальное расстояние между частицами в x- и y-направлениях равно L/2, а не L. Численный алгоритм. Краевые условия

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

Численный алгоритм. Краевые условия Рис. 1. Пример периодических краевых условий в двумерном случае. Обратите внимание на то, что частица 1 собирается покинуть левую грань центральной клетки и войти в центральную клетку справа. Правило ближайшей частицы означает, что расстояние между частицами 1 и 2 определяется жирной линией.

Это правило означает, что частица 1 из центральной клетки взаимодействует только с ближайшим экземпляром частицы 2; взаимодействие полагается равным нулю, если расстояние до копии больше L/2 (см. рис.1). Поскольку из данного правила следует, что мы можем визуально представлять себе центральную клетку в виде тора, это использование периодических краевых условий вместе с правилом ближайшей частицы было бы точнее назвать тороидальными краевыми условиями. Однако мы верны принятым традициям и называем эти условия периодическими краевыми условиями. Отметим, что использование периодических краевых условий означает, что все узлы ящика эквивалентны. Численный алгоритм. Краевые условия

Лекция окончена Нажмите клавишу для выхода