Моделирование микроканонического ансамбля методом Монте Карло в пакете MATLAB (стр. 2 )

Задача 1.
Проведите численные расчеты для заданных значений числа частиц (N=10, 20, 40, 80), и различных значений

отношения

. Изобразите графически зависимости коэффициента принятия (Accept) от

при заданном числе частиц. Используя метод наименьших квадратов, подберите аналитическое

выражение, описывающее зависимость коэффициента принятия (Accept) от

. При каком отношении

значение коэффициента принятия оказывается 0.5?

Задача 2.
Модифицируйте функцию Demon, представленную на рис. 2, так чтобы номера частиц, скорость которых подлежит изменению (переменная Ip), выбирались последовательно. Приводит ли изменение способа выбора номера частицы к изменению средней энергии системы на частицу и средней энергии демона?
Задача 3.
Можно ввести кинетическое определение температуры, в соответствие с которым, температура определяется из соотношения

(3)

где

- средняя кинетическая энергия на частицу, Tkin температура газа, измеряемая в единицах kB.

Модифицируйте функцию Demon( ), представленную на рис. 2, так чтобы в программе дополнительно вычислялось и дополнительно выводилось в векторе, возвращаемом данной функцией, значение средней кинетической энергии на одну частицу. Используя соотношение (3), получите Tkin. Как связаны значения температуры и энергии демона?
Один из наиболее важных результатов, полученных в задаче 3, состоит в демонстрации возможности введения кинетического определения температуры, как средней на одну частицу кинетической энергии частицы микроканонического ансамбля. При этом значение температуры с некоторой погрешностью оказывается равной энергии демона. В то же время микроканонический ансамбль нельзя считать абсолютно адекватной моделью реальных статистических систем, так данные системы не являются термодинамически изолированными, но находятся в тепловом контакте с окружающей средой. Это приводит к возникновению потока тепла между изучаемой системой и окружающей средой. Так как размер исследуемой системы обычно много меньше размера окружающей среды, первую будем называть микросистемой, а последнюю тепловым резервуаром (термостатом). При этом закон сохранения полной энергии относится к составной системе, состоящей из микросистемы и теплового резервуара.
Рассмотрим большое число воображаемых копий микросистемы и теплового резервуара. Рассматриваемые как единое целое микросистема и тепловой резервуар являются изолированной термодинамической системой и могут быть описаны с помощью микроканонического ансамбля. Так как в данном случае наибольший интерес представляют равновесные значения физических величин, описывающих микросистему, необходимо знать вероятность Ps, с которой микросистема обнаруживается в состоянии s с энергией Es. Ансамбль, который описывает распределение вероятностей состояний микросистемы, находящейся в термодинамическом равновесии с тепловым резервуаром, называется каноническим.
В общем случае в качестве микросистемы, может выступать любая макроскопическая система, размер которой меньше теплового резервуара, в том числе и отдельная частица, выделенная среди других частиц теплового резервуара. Примером такой микросистемы служит демон, которого можно считать микросистемой, микроканоническое состояние которой определяется только ее энергией. Таким образом, метод отыскания распределения вероятностей в каноническом ансамбле заключается в том, чтобы выполнить численное моделирование демона, который обменивается энергией с идеальным газом, состоящим из N частиц. Идеальный газ в данном случае играет роль термостата и задача состоит в определении плотности вероятности P(Ed) того, что демон имеет энергию Ed. Вычислительный алгоритм, позволяющий решить данную задачу, достаточно очевиден:

Провести моделирование канонического ансамбля методом Монте-Карло, сохраняя на каждом шаге Монте-Карло значение энергии демона. Для найденной реализации значений энергии демона вычислить распределение вероятностей P(Ed)E. Определить функциональную зависимость P(Ed)E от температуры.

Для реализации данного алгоритма в пакете MATLAB создадим m-файл DemonE.m, содержащий описание функции, возвращающей мгновенные (на каждом шаге Монте-Карло) значения энергии демона.

% листинг файла DemonE. m
function z = DemonE(N,Esystem,NTrial,Vel,dV)
% функция, возвращающая мгновенные значения энергии демона
% N число частиц системы
% Esystem энергия системы
% NTrial число испытаний
% Vel вектор, возвращенный функцией InitD
% dV максимальное изменение скорости
Edemon=0; z(1)=Edemon; k=1; for j=1:NTrial  for i=1:N  dv=(2*rand(1)-1)*dV; % случайное изменение скорости
 Ip=floor(N*rand(1)+1); % случайный выбор номера частицы
 Vtrial=Vel(Ip)+dv;  de=0.5*(Vtrial.^2-Vel(Ip).^2); % пробное изменение энергии
 if de<=Edemon % если энергия уменьшается, то пробное изменение
 % принимается
 Vel(Ip)=Vtrial;  Edemon=Edemon-de;  end;  k=k+1;  z(k)=Edemon;  end; end;

Далее необходимо выполнить следующую последовательность команд:

>> Esystem=40; % энергия системы
>> N=40; % число частиц системы
>> NTrial=500; % число испытаний
>> dV=2*2^0.5; % максимальное изменение скорости
>> Vel=InitD(N,Esystem); % вычисление начальных скоростей частиц
>> Ed=DemonE(N,Esystem,NTrial,Vel,dV);% вычисление мгновенных
% значений энергии демона
% визуализация зависимости мгновенных значений энергии демона
% от времени
>> i=1:length(Ed);
>> figure(1);plot(i,Ed,'k');axis([0 20000 0 8])
% вычисление и визуализация распределения вероятностей
>> x1=min(Ed);
>> x2=max(Ed);
>> Nint=50;
>> i=1:Nint;
>> x(i)=x1+(x2-x1)/(Nint-1)*(i-1);
>> h=hist(E,x);
>> figure(2);bar(x,h); colormap white % вычисление и визуализация функции, аппроксимирующей
% распределение
>> F10=inline('u(1)*exp(-u(2)*z)','u','z') % задание аппроксимирующей функции
F10 =  Inline function:  F10(u,z) = u(1)*exp(-u(2)*z)
>> beta=nlinfit(x,h,F10,[3000 1]); % вычисление коэффициентов аппроксимирующей функции
% визуализация гистограммы последовательности мгновенных значений
% энергии демона и аппроксимирующей функции
>> Ni=500;
>> i=1:Ni;
>> X(i)=x1+(x2-x1)/(Ni-1)*(i-1);
>> figure(3);bar(x,h); colormap white
>> hold on
>> plot(X(i),F10(beta,X(i)),'k')
>> hold off