Всем привет 👋 В прошлой части этой серии статей я закончил на том, что в нашей симуляции не выполняется закон сохранения энергиисохраняется энергия из-за какой-то проблемы. Давайте разбираться, что же не так.

Думаю, что в модели проблемы нет, ведь она известна уже долгое время и часто используется в подобных расчётах. Тогда получается, что проблема в нашей программе. Сейчас алгоритм работает следующим образом:

Задаётся некоторый временной шаг (например, 1 мс). Далее в цикле выполняются следующие действия:

  1. Рассчитывается ускорение, исходя из силы (это точно верно, ведь в потенциале Леннарда-Джонса и втором законе Ньютона мы не сомневаемся)
  2. Ускорение умножается на время для получения изменения скорости
  3. Скорость умножается на время для получения изменения координаты

Если ошибка не в первом шаге, то она кроется во 2 или 3. Давайте прочитаем, каким образом связано ускорение и координата тела. На википедии мы находим следующее:

Pasted Graphic.tiff

Это немного отличается от того, что мы используем в нашей симуляции:

Pasted Graphic 1.tiff

Возможно в этом и кроется проблема! Если изучить вторую формулу, то можно понять, что она верна лишь для того случая, когда ускорение всё время t равняется a, то есть в случае равноускоренного движения. Действительно,

Pasted Graphic 4.tiff

если ускорение все время постоянно (точки над переменной в физике обозначают производную; 1 точка - первая производная, 2 точки - вторая производная…).

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

Что же нам поможет? Давайте сначала разберёмся с математикой, а потом с тем, как её реализовать в нашей программе.

🧮 Математика

Например, для двух частиц у нас получится следующая система:

Pasted Graphic_1.tiff

Где F - функция, возвращающая нам силу. То есть получается довольно интересная система уравнений, решением которой являются не числа, как это бывает в обычных уравнениях, а функции! Получается, что для r1 и r2 мы должны подобрать такие функции, что их вторые производные будут удовлетворять данным уравнениям.

Такие уравнения называются дифференциальными, и в общем случае решение этих уравнений - это невероятно сложная задача. Некоторые важнейшие дифференциальные уравнения до сих пор остаются без решений, а для каких-то из них доказано, что в общем случае решений нет (например, https://ru.wikipedia.org/wiki/Задача_трёх_тел).

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

Например, есть следующее уравнение:

Pasted Graphic 3.tiff

Его решением является функция:

Pasted Graphic 4_1.tiff

Проверить, что это является решением очень просто: необходимо просто подставить нашу функцию в уравнение. Давайте попробуем это сделать:

Pasted Graphic 5.tiff

Но можно заметить, что здесь появляется несколько новых переменных (A и φ0). Так происходит часто, иногда решение определяется с точностью до константного множителя, иногда - с точностью до константы. Иногда, чтобы убрать эти переменные, задают начальные значения. Например, для этой задачи мы можем задать значение x(0).

Вернёмся теперь к нашей изначальной системе уравнений. Она сложнее предыдущей задачи, ведь у нас 2 связанных между собой уравнения. И это только для двух частиц! Представим, что мы хотим смоделировать 1000 частиц, у нас получится 1000 дифференциальных уравнений, каждое из которых будет связано с остальными. Здесь уже не придумаешь явную функцию, отвечающую всем этим условиям, как в предыдущей задаче. Но, вспоминая что у нас есть компьютер, мы задумываемся, а не можем ли мы переложить на него сложную задачу по решению этих уравнений. И ответ - можем!

🧑‍💻 Программирование

Для решения нашей задачи существует область математики под названием «Вычислительная математика». Она изучает методы вычисления различных математических выражений. Некоторые такие методы позволяют находить решения для систем дифференциальных уравнений. Но, к сожалению, многие из них не дают нам точного решения, то есть не позволяют найти функции, которые бы удовлетворяли всем уравнениям системы. Однако, несмотря на этот огромный недостаток, есть и огромный плюс: мы можем находить приблизительные решения для огромного класса дифференциальных уравнений. Более того, есть методы, которые позволяют найти решение с любой необходимой нам точностью, повторяя одни и те же операции много раз. Это то, что нам и нужно, ведь компьютер очень хорошо и быстро выполняет строго заданные ему операции. Из всего многообразия алгоритмов я решил выбрать метод Рунге-Кутты. Это один из тех методов, который позволяет достигать любой точности. Для моих задач подойдёт метод Рунге-Кутты четвёртого порядка (можно взять больший порядок, чтобы увеличить точность, но и длительность вычислений тоже возрастёт).

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

Pasted Graphic 6.tiff

То для вычисления следующих значений нам важны лишь значения функций f в конкретных точках, то есть функции f могут быть любыми!

Остаётся последняя проблема - у нас-то в задаче уравнения второго порядка (в них присутствует вторая производная), а здесь первого. Но это не проблема, ведь любое дифференциальное уравнение второго порядка можно свести к системе дифференциальных уравнений первого порядка:

Pasted Graphic 8.tiff

Таким образом, для системы из 1000 частиц мы получаем систему из 2000 дифференциальных уравнений первого порядка. Это тяжело даже записать человеку, но компьютер может производить вычисления с подобными уравнениями десятки раз в секунду.

Но я так и не показал, как метод Рунге-Кутты (4 порядка) работает. Вот основные формулы:

Если у нас дана такая задача:

Pasted Graphic_2.tiff

То значение в следующей точке мы можем найти следующим образом:

Pasted Graphic 1_1.tiff

Здесь h - это тот самый шаг, который мы можем сами выбирать.

Получается, что в нашем случае такие операции нам придётся проделывать для каждого уравнения в системе.

Исходя из всего этого, я переделал алгоритм из предыдущей статьи. Для начала посмотрим, что будет происходить с системой из двух частиц:

ezgif-3-a295316202.gif

Визуально результат очень похож на симуляцию из предыдущей статьи, но решилась ли проблема с энергией? И ответ - да. Если в предыдущей статье отклонения от теоретической энергии составляли где-то 3% от теоретического значения, то здесь отклонение составляет 0.000000001%! Это во много раз лучше, чем у нас было до этого. 

Давайте посмотрим на систему из большего числа частиц.

Particles.gif

Если в предыдущей версии программы изменения энергии для этой системы составляли более 10% от абсолютного значения, то теперь составляют около десятой процента.

Для наглядного сравнения приведу график зависимости энергии в системе из 25 частиц от времени. Здесь синий цвет - это теоретическое значение, зелёный - значение, полученное при расчётах с дифференциальными уравнениями, а серый - без них. На оси X отложено количество итераций, на оси Y - значение энергии в условных единицах. (количество итераций, начальные условия и шаг были одинаковыми при обоих запусках)

Снимок экрана 2022-09-09 в 21.30.16.png

При изучении этого графика может возникнуть вопрос: почему в первые сотни итераций энергия в симуляциях очень сильно отклоняется от теоретического значения, а затем выравнивается и идёт параллельно с теоретическим значением? Думаю, что ответ кроется в том, что наибольшая ошибка в симуляциях возникает, когда функция силы от расстояния между частицами меняется очень резко, а это происходит, когда атомы находятся очень близко друг к другу. Так как я изначально расположил атомы случайно, то какие-то из них оказались на таком расстоянии, что сила менялась очень быстро при их удалении. Поэтому обычный метод выдал очень большое отклонение в самом начале, а после этого, когда атомы разлетелись на достаточные расстояния, где значения меняются не так резко, энергия практически перестала меняться. Программа с дифференциальными уравнениями гораздо лучше справляется в таких крайних условиях, поэтому и ошибка получилось меньше.

Добиться ещё большей точности в симуляции мы можем с помощью уменьшения шага, но это, как я много раз повторял, ведёт к увеличению времени вычислений, что не всегда так хорошо. Таким образом, приходится балансировать между точностью и временем вычислений, чтобы получить что-то, что возможно использовать. Для наглядности вот график зависимости медианы модуля ошибки в зависимости от шага. На оси Y здесь отложена медиана модуля ошибки в логарифмическом масштабе. По оси X отложен шаг

Снимок экрана 2022-09-10 в 00.50.51.png

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

🤔 Заключение

Теперь у нас есть программа, позволяющая моделировать различные системы частиц, а значит теперь мы можем с ними экспериментировать прямо на компьютере! Конечно, на идеальную точность претендовать не получится, но посмотреть, что выйдет из этих симуляций всё равно интересно. В следующей части этой серии статей мы проверим законы молекулярно-кинетической теории, а если всем будет интересно, то можем и до аэродинамики коровы дойти. Если кому-то интересно разобраться в том, что за программу я написал, то вот ссылка на репозиторий с проектом: https://github.com/MishkaSimakov/LennardJonesPotentialSimulator. Буду рад, если вы найдёте какие-то ошибки в нём и напишите мне об этом. Спасибо за чтение!

📚 Литература



Нет комментариев