Моделирование межмолекулярного взаимодействия (часть 2)
Всем привет 👋 В прошлой части этой серии статей я закончил на том, что в нашей симуляции не выполняется закон сохранения энергиисохраняется энергия из-за какой-то проблемы. Давайте разбираться, что же не так.
Думаю, что в модели проблемы нет, ведь она известна уже долгое время и часто используется в подобных расчётах. Тогда получается, что проблема в нашей программе. Сейчас алгоритм работает следующим образом:
Задаётся некоторый временной шаг (например, 1 мс). Далее в цикле выполняются следующие действия:
- Рассчитывается ускорение, исходя из силы (это точно верно, ведь в потенциале Леннарда-Джонса и втором законе Ньютона мы не сомневаемся)
- Ускорение умножается на время для получения изменения скорости
- Скорость умножается на время для получения изменения координаты
Если ошибка не в первом шаге, то она кроется во 2 или 3. Давайте прочитаем, каким образом связано ускорение и координата тела. На википедии мы находим следующее:
Это немного отличается от того, что мы используем в нашей симуляции:
Возможно в этом и кроется проблема! Если изучить вторую формулу, то можно понять, что она верна лишь для того случая, когда ускорение всё время t равняется a, то есть в случае равноускоренного движения. Действительно,
если ускорение все время постоянно (точки над переменной в физике обозначают производную; 1 точка - первая производная, 2 точки - вторая производная…).
Но ведь в нашем случае, какое бы мы маленькое t не взяли, нельзя считать, что всё это время ускорение постоянно, ведь наши частицы постоянно меняют свои положения, а значит и силы, действующие на них, меняются постоянно, каждое мгновение.
Что же нам поможет? Давайте сначала разберёмся с математикой, а потом с тем, как её реализовать в нашей программе.
🧮 Математика
Например, для двух частиц у нас получится следующая система:
Где F - функция, возвращающая нам силу. То есть получается довольно интересная система уравнений, решением которой являются не числа, как это бывает в обычных уравнениях, а функции! Получается, что для r1 и r2 мы должны подобрать такие функции, что их вторые производные будут удовлетворять данным уравнениям.
Такие уравнения называются дифференциальными, и в общем случае решение этих уравнений - это невероятно сложная задача. Некоторые важнейшие дифференциальные уравнения до сих пор остаются без решений, а для каких-то из них доказано, что в общем случае решений нет (например, https://ru.wikipedia.org/wiki/Задача_трёх_тел).
Но перед тем, как подступаться к нашей системе, давайте посмотрим на какое-нибудь простое дифференциальное уравнение и его решение, чтобы подробнее в них разобраться.
Например, есть следующее уравнение:
Его решением является функция:
Проверить, что это является решением очень просто: необходимо просто подставить нашу функцию в уравнение. Давайте попробуем это сделать:
Но можно заметить, что здесь появляется несколько новых переменных (A и φ0). Так происходит часто, иногда решение определяется с точностью до константного множителя, иногда - с точностью до константы. Иногда, чтобы убрать эти переменные, задают начальные значения. Например, для этой задачи мы можем задать значение x(0).
Вернёмся теперь к нашей изначальной системе уравнений. Она сложнее предыдущей задачи, ведь у нас 2 связанных между собой уравнения. И это только для двух частиц! Представим, что мы хотим смоделировать 1000 частиц, у нас получится 1000 дифференциальных уравнений, каждое из которых будет связано с остальными. Здесь уже не придумаешь явную функцию, отвечающую всем этим условиям, как в предыдущей задаче. Но, вспоминая что у нас есть компьютер, мы задумываемся, а не можем ли мы переложить на него сложную задачу по решению этих уравнений. И ответ - можем!
🧑💻 Программирование
Для решения нашей задачи существует область математики под названием «Вычислительная математика». Она изучает методы вычисления различных математических выражений. Некоторые такие методы позволяют находить решения для систем дифференциальных уравнений. Но, к сожалению, многие из них не дают нам точного решения, то есть не позволяют найти функции, которые бы удовлетворяли всем уравнениям системы. Однако, несмотря на этот огромный недостаток, есть и огромный плюс: мы можем находить приблизительные решения для огромного класса дифференциальных уравнений. Более того, есть методы, которые позволяют найти решение с любой необходимой нам точностью, повторяя одни и те же операции много раз. Это то, что нам и нужно, ведь компьютер очень хорошо и быстро выполняет строго заданные ему операции. Из всего многообразия алгоритмов я решил выбрать метод Рунге-Кутты. Это один из тех методов, который позволяет достигать любой точности. Для моих задач подойдёт метод Рунге-Кутты четвёртого порядка (можно взять больший порядок, чтобы увеличить точность, но и длительность вычислений тоже возрастёт).
Поговорим немного об этом методе. Он позволяет находить приблизительные решения для систем дифференциальных уравнений первого порядка. При этом следующие значения считаются, исходя из предыдущих. Мы задаём такой параметр, как шаг. Чем он меньше, тем больше точность, но и больше времени потребуется для вычисления одного и того же временного промежутка. Что самое интересное, если у нас есть такая система:
То для вычисления следующих значений нам важны лишь значения функций f в конкретных точках, то есть функции f могут быть любыми!
Остаётся последняя проблема - у нас-то в задаче уравнения второго порядка (в них присутствует вторая производная), а здесь первого. Но это не проблема, ведь любое дифференциальное уравнение второго порядка можно свести к системе дифференциальных уравнений первого порядка:
Таким образом, для системы из 1000 частиц мы получаем систему из 2000 дифференциальных уравнений первого порядка. Это тяжело даже записать человеку, но компьютер может производить вычисления с подобными уравнениями десятки раз в секунду.
Но я так и не показал, как метод Рунге-Кутты (4 порядка) работает. Вот основные формулы:
Если у нас дана такая задача:
То значение в следующей точке мы можем найти следующим образом:
Здесь h - это тот самый шаг, который мы можем сами выбирать.
Получается, что в нашем случае такие операции нам придётся проделывать для каждого уравнения в системе.
Исходя из всего этого, я переделал алгоритм из предыдущей статьи. Для начала посмотрим, что будет происходить с системой из двух частиц:
Визуально результат очень похож на симуляцию из предыдущей статьи, но решилась ли проблема с энергией? И ответ - да. Если в предыдущей статье отклонения от теоретической энергии составляли где-то 3% от теоретического значения, то здесь отклонение составляет 0.000000001%! Это во много раз лучше, чем у нас было до этого.
Давайте посмотрим на систему из большего числа частиц.
Если в предыдущей версии программы изменения энергии для этой системы составляли более 10% от абсолютного значения, то теперь составляют около десятой процента.
Для наглядного сравнения приведу график зависимости энергии в системе из 25 частиц от времени. Здесь синий цвет - это теоретическое значение, зелёный - значение, полученное при расчётах с дифференциальными уравнениями, а серый - без них. На оси X отложено количество итераций, на оси Y - значение энергии в условных единицах. (количество итераций, начальные условия и шаг были одинаковыми при обоих запусках)
При изучении этого графика может возникнуть вопрос: почему в первые сотни итераций энергия в симуляциях очень сильно отклоняется от теоретического значения, а затем выравнивается и идёт параллельно с теоретическим значением? Думаю, что ответ кроется в том, что наибольшая ошибка в симуляциях возникает, когда функция силы от расстояния между частицами меняется очень резко, а это происходит, когда атомы находятся очень близко друг к другу. Так как я изначально расположил атомы случайно, то какие-то из них оказались на таком расстоянии, что сила менялась очень быстро при их удалении. Поэтому обычный метод выдал очень большое отклонение в самом начале, а после этого, когда атомы разлетелись на достаточные расстояния, где значения меняются не так резко, энергия практически перестала меняться. Программа с дифференциальными уравнениями гораздо лучше справляется в таких крайних условиях, поэтому и ошибка получилось меньше.
Добиться ещё большей точности в симуляции мы можем с помощью уменьшения шага, но это, как я много раз повторял, ведёт к увеличению времени вычислений, что не всегда так хорошо. Таким образом, приходится балансировать между точностью и временем вычислений, чтобы получить что-то, что возможно использовать. Для наглядности вот график зависимости медианы модуля ошибки в зависимости от шага. На оси Y здесь отложена медиана модуля ошибки в логарифмическом масштабе. По оси X отложен шаг
Видно, что в среднем уменьшение шага приводит к уменьшению ошибки. Огромные значения ошибки при шаге = 1 связаны с тем, что это слишком большой шаг, не подходящий даже для самых простых симуляций, ведь при таких больших значениях частицы начинают пролетать сквозь друг друга, разлетаться с слишком большими скоростями и т.п.
🤔 Заключение
Теперь у нас есть программа, позволяющая моделировать различные системы частиц, а значит теперь мы можем с ними экспериментировать прямо на компьютере! Конечно, на идеальную точность претендовать не получится, но посмотреть, что выйдет из этих симуляций всё равно интересно. В следующей части этой серии статей мы проверим законы молекулярно-кинетической теории, а если всем будет интересно, то можем и до аэродинамики коровы дойти. Если кому-то интересно разобраться в том, что за программу я написал, то вот ссылка на репозиторий с проектом: https://github.com/MishkaSimakov/LennardJonesPotentialSimulator. Буду рад, если вы найдёте какие-то ошибки в нём и напишите мне об этом. Спасибо за чтение!
📚 Литература
- https://www.mccme.ru/free-books/dubna/anosov.pdf (интересная книжка для школьников про дифференциальные уравнения)
- http://mathprofi.ru/metody_eilera_i_runge_kutty.html (про методы приближённого решения дифференциальных уравнений)