Метод Монте-Карло

Метод Монте-Карло (МК) был разработан в конце второй мировой войны Дж. фон Нойманом, С. Уламом и Н. Метрополисом для изучения диффузии нейтронов в делящихся материалах. Датой рождения метода МК считается дата появления статьи с таким названием.[1] Название метода связано с тем, что в нем широко используются случайные числа (Монте-Карло, как известно, славится своими игорными заведениями, где случайность – основа всех игр, притягивающая азартных людей в это княжество). Моделирование основано на генерации случайных чисел с помощью арифметических и логических операций, которых к настоящему времени разработано большое множество. На основе этих чисел с помощью компьютера разыгрываются случайные реализации тех или иных процессов, состояний и событий. Метод МК может быть использован для моделирования не только случайных событий, но также для вполне детерминистских расчетов. Например, с его помощью можно провести численное интегрирование сложных функций. В качестве примера рассмотрим расчет площади сектора круга и определение числа p.

7.1. Расчет числа p с помощью метода Монте Карло

Рассмотрим круг единичного радиуса. Число p может быть рассчитано как площадь этого круга. Последняя равна учетверенной площади закрашенной области на рис.1 – четверти круга. С другой стороны, площадь квадрата ABCO равна единице. Число p равно

, (7.1)

На квадрате ABCO генерируются равномерно-случайно распределенные точки. То есть, для каждой точки генерируются две случайные координаты, являющиеся реализациями случайной величины, равномерно распределенной в отрезке [0,1). Определяется расстояние каждой точки до начала координат. Если это расстояние меньше 1, то точка принадлежит к четверти круга; из полного числа розыгрышей определяется число попаданий в четверть круга . Ясно, что отношение будет равно отношению площадей четверти круга и квадрата, то есть площади четверти круга. Следовательно, число p может быть рассчитано как

. (7.2)

Разыграв большое число случайных точек, можно с весьма большой точностью определить значение числа p. Ключевым моментом в этом методе является генерация случайных чисел из равномерного распределения. Это является основой метода МК.

 

Рис. 6. К расчету числа p методом Монте Карло

Метод МК имеет преимущества перед обычными, детерминистскими методами при расчете интегралов в многомерном пространстве. Вычисление таких интегралов необходимо в статистической физике для расчета термодинамических величин систем многих частиц, которые определяются как средние некоторых величин по ансамблю. Именно для расчета таких средних и был разработан алгоритм метода МК, который носит название алгоритма Метрополиса.

7.2. Средние значения в статистической физике

Рассмотрим систему, состоящую из N частиц, описываемую функцией Гамильтона , где ‑ радиусы-векторы и импульсы частиц. Состояние системы определяется точкой фазового пространства, осями которого являются величин – координаты и проекции импульсов всех частиц. Множество всех возможных состояний системы определяет доступную ей область фазового пространства, которую обозначим W. Как известно, среднее значение величины A, являющейся функцией состояния системы, определяется интегралом

, (7.3)

где Z ‑ статистический вес системы:

, (7.4)

Рассмотрим случай, когда функция А зависит только от конфигурации системы, но не от импульсов частиц: . Гамильтониан системы состоит из суммы кинетической и потенциальной энергий:

. (7.5)

Интегралы (7.3) и (7.4) могут быть разделены на интегралы по пространственным координатам и импульсам, а для А, зависящих только от координат, и в числителе, и в знаменателе получим один и тот же интеграл

. (7.6)

В этом случае расчет средней величины сводится к интегрированию по пространственным координатам:

, (7.7)

где ‑ конфигурационный интеграл, равный

. (7.8)

Величина

(7.9)

определяет плотность вероятности нахождения системы в конфигурации . Поэтому (7.7) можно переписать в виде

. (7.10)

Можно попытаться посчитать интеграл (7.10) с использованием квадратурных методов, например, метода Симпсона. Для этого отрезок интегрирования по каждой координате необходимо поделить на некоторое достаточно большое число m маленьких отрезков. Следовательно, для расчета интеграла необходимо посчитать значения интегрируемой функции в точках. При расчете 10-кратного интеграла с делением на m =100 точек это число равно , что уже представляет собой астрономическое число. Отсюда видим, что обычные методы интегрирования для расчета средних величин при большом числе частиц N неприменимы.

Можно также попытаться рассчитать интеграл методом Монте Карло, как в примере расчета числа p. Для этого нужно сделать следующее:

1) Генерируя 3 N случайных координат, найти одно состояние системы;

2) Рассчитать больцмановский множитель и произведение ;

3) Прибавить рассчитанные величины к накопленным предыдущими шагами суммам и вернуться к п.1.

4) После определенного количества итераций N _ит определить среднее значение как

. (7.11)

Этот подход, однако, также невыполним ввиду его большой трудоемкости. В суммы выражения (7.11) будет входить огромное множество состояний с пренебрежимо малым множителем Больцмана, которым соответствует высокая энергия. В суммы главный вклад вносят состояния только небольшой области конфигурационного пространства, в которых указанный множитель имеет достаточно большое значение. По трудоемкости этот метод не будет отличаться существенно от метода квадратур.

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

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

7.3. Алгоритм Метрополиса

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

Пусть система находится в конфигурации , которую мы обозначим о (old), имеющей ненулевую вероятность . Из этого состояния путем случайных изменений координат генерируется новое состояние n (new). Вероятность нового состояния равна . Для построения алгоритма Метрополиса необходимо определить вероятность перехода из старого состояния в новое , то есть нужно разработать правило для принятия решения, допускать переход или отказаться.

Чтобы понять это правило, представим себе, что одновременно проводится очень большое число М моделирований Монте Карло, где М намного больше, чем полное число возможных конфигураций системы. Число систем в любой конфигурации о обозначим m (o); это число пропорционально вероятности состояния . Нас интересует равновесное распределение состояний и переходов между ними. При таком равновесии среднее число переходов , выводящих систему из состояния о, равно числу переходов в это состояние со всех других. Можно наложить и более сильное условие, что в равновесии среднее число принятых переходов из состояния о в любое другое состояние n равно среднему числу обратных переходов из состояния n в о. Таким образом, цепочка состояний, которую нужно генерировать, должна удовлетворять следующему условию детального равновесия

. (7.12)

Этому условию можно удовлетворить различными способами. В алгоритме Метрополиса принимается следующее правило:

(7.13)

Видно, что это правило удовлетворяет условию (7.12).

Правило (7.13) позволяет генерировать, начиная с некоторого начального состояния, требуемое количество m состояний, в соответствии с их вероятностью, и по ним произвести усреднение физических величин:

. (7.14)

Поскольку состояния генерируются в соответствии с их вероятностью, умножения на фактор Больцмана в этом выражении нет.

Алгоритм Метрополиса может быть реализован следующим образом.

В системе случайным образом выбирается одна частица. Пусть ее текущие координаты равны . Этой частице придается случайное перемещение в случайном направлении:

,

, (7.15)

,

где ‑ максимальное отклонение в каждом направлении. Для изменения каждой координаты генерируется новое случайное число x. Для новой конфигурации рассчитывается энергия системы; для этого достаточно пересчитать энергии атомов небольшой области вблизи выбранной, так как потенциал взаимодействия короткодействующий,‑ это можно сделать введением списка соседей, как и в МД. Энергия новой конфигурации сравнивается с энергией старой . Если , новая конфигурация принимается и становится стартовой для следующего шага. Если же , то рассчитывается отношение факторов Больцмана нового и старого состояний и генерируется случайное число rand(0,1), равномерно распределенное на отрезке (0,1). Если

, (7.16)

новая конфигурация принимается; в противном случае она отвергается, и старая конфигурация остается стартовой для новой попытки.

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

Атом, смещение которого моделируется, может быть выбран на каждом шаге случайным образом; можно смещать атомы последовательно, друг за другом. На каждом шаге могут быть смещены несколько атомов. Путем оптимизации параметра в этом случае можно более эффективно исследовать конфигурационное пространство.

Так же, как при МД моделировании, моделирование методом МК требует определенного числа шагов для достижения равновесия. Независимо от исходной конфигурации, алгоритм Метрополиса через некоторое, достаточно большое число шагов приводит к равновесию, поскольку начальное состояние в Марковской цепи «забывается». В течение этих шагов можно отслеживать изменения энергии, структурных параметров, среднеквадратичных отклонений и т.д. Средние значения определяются с использованием шагов после достижения равновесия.

Так же, как в МД, в методе МК моделирование можно проводить для различных ансамблей (микроканонических, канонических, изотермо-изобарических, больших канонических). Главным отличием метода МК от МД является отсутствие времени и импульсов, поэтому с его помощью могут быть определены только средние значения величин, зависящих только от координат атомов.

7.4. Моделируемый отжиг

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

Процедура моделируемого отжига основана на рассмотренном выше алгоритме Метрополиса. Рассмотрим систему из N частиц, потенциальную энергию которой надо минимизировать. Для этого выбирается некоторая исходная конфигурация. Начиная с этой конфигурации, генерируются новые состояния по алгоритму Метрополиса при некоторой достаточно высокой температуре T. Далее температура постепенно уменьшается до указанной температуры (можно до очень низкой, близкой к нулю температуры). При высокой температуре алгоритм Метрополиса допускает генерацию значительного количества состояний с большим увеличением энергии, поэтому будут исследованы области конфигурационного пространства, весьма далекие от положения минимума энергии. С уменьшением температуры таких изменений конфигурации будет все меньше, и поиск все больше сосредоточивается вблизи минимума. Последние 40% шагов делается вблизи минимума. Это типично для алгоритма моделируемого отжига.

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

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

7.5. Генераторы случайных чисел

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

Как правило, генераторы моделируют случайные числа, имеющие равномерное распределение в интервале (0,1). Если генерировать точки в N -мерном пространстве, для каждой координаты точки необходимо одно случайное число, то есть, для каждой точки нужно генерировать N случайных чисел. Большое число точек, генерированных в N -мерном кубе, должно быть равномерно распределено в объеме этого куба. Не любой алгоритм генерации удовлетворяет этому минимальному требованию. Часто точки имеют тенденцию лежать на гиперповерхности в этом пространстве.

Большинство генераторов случайных чисел основано на линейном конгруэнтном методе. Каждое текущее число определяется из предыдущего путем умножения его на некоторое постоянное число a, сложения результата с другим постоянным числом b и нахождения остатка от деления на третье постоянное число m. Вся последовательность определяется однозначно начальным числом . Таким образом:

, (7.17)

Этот алгоритм генерирует целые числа в интервале от 1 до m‑ 1, поделив которые на m, можно получить вещественные числа в интервале (0,1). Генерируемая последовательность будет периодической с периодом, равным m, поэтому число m должно быть достаточно большим.

Дата добавления: 2015-10-02