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

(10)

Предваряя применения алгоритма демона к исследованию модели Изинга методом микроканонического ансамбля, необходимо получить выражение, связывающее энергию демона и температуру термостата. Напомним, что в непрерывном случае распределение энергии демона подчинялось формуле Больцмана (6). Предположим, что данное распределение вероятностей справедливо для любой макроскопической системы, находящейся в состоянии термодинамического равновесия. Тогда <Ed> равно

(11)

где суммы вычисляются по всем возможным значениям Ed. Минимальная ненулевая потеря энергии системы в нулевом магнитном поле составляет, как видно из рис. 5, 2sJ, где s - суммарный спин ближайших соседей, опрокидывающегося спина.

Рис. 5. Изменение энергии, обусловленное изменения ориентации центрального спина

В одномерном случаем суммарный спин ближайших соседей равен 0, или 2, т. е. минимальная не нулевая потеря энергии равняется 2J. Следовательно, энергия демона, может равняться 0, 2J, 4J, ... . Если ввести обозначение x=2J/T, то выражение (11) для бесконечной решетки принимает вид

(12)

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

	(13)
	(14)

Разделив (13) на (14) и выполнив очевидное преобразование, получаем

(15)

Подставив в (15) вместо x выражение 2J/T и решив получившее уравнение относительно T, получаем

(16)

Задача 5.
Получите аналитические выражения для средней энергии демона для системы, состоящей из конечного числа спинов. Может ли быть найдено аналитическое выражение для зависимости температуры от энергии демона в данном случае? Сравните значения средней энергии демона для системы, состоящей из N спинов <Ed>N, получаемые по точной формуле с соответствующими значениями <Ed>, получаемыми по формуле (15). Постройте зависимости |<Ed>N-<Ed>| при различных температурах. Оцените погрешность определения значения температуры, получаемого по формуле (16), для системы, состоящей из N спинов. Как зависит данная погрешность от N и T?

Отметим, что в одномерной модели Изинга демон должен выбирать спины случайно, что позволит избежать подсчета периодически повторяющихся конфигураций. Так как нашей целью является оценка термодинамических характеристик бесконечной термодинамической системы, следует учесть краевые условия. В качестве таковых в описываемом ниже документе нами выбраны периодические (тороидальные) краевые условия: решетка считается кольцом, в котором спины si, находящиеся в узлах i=1 и i=N, взаимодействуют друг с другом. Это обеспечивает равное число взаимодействий для всех спинов в изучаемой системе.
Для моделирования одномерной модели Изинга методом микроканонического ансамбля в пакете MATLAB необходимо создать файл Ising.m, содержащий описание функции, возвращающей значения полной энергии системы, энергии демона, намагниченности и среднее число принятия решений.

% листинг файла Ising. m
function [Es,Ed,SpM,Accept]=Ising(Nspin,J,h,Esi,NTrial)
% функция, возвращающая мгновенные значения полной энергии системы
% (Es), энергии демона (Ed), намагниченность (SpM) и среднее число
% принятия (Accept)
% Nspin число спинов системы
% J константа обменного взаимодействия
% h внешнее магнитное поле
% Esi конечная энергия системы
% NTrial число испытаний
% задание конфигурации спинов в момент времени t = 0
for i=1:Nspin  s(i)=1; end; M=Nspin; Esystem=-(J+h)*Nspin; Edemon=2*J*ceil((Esi-Esystem)/(2*J)); Es(1)=Esystem; % энергия системы в момент времени t = 0
Ed(1)=Edemon; % энергия демона в момент времени t = 0
SpM(1)=M; % магнитный момент системы в момент времени t = 0
Accept=0; k=1;
% реализация метода микроканонического ансамбля
for i=1:NTrial  for j=1:Nspin  Ispin=floor(Nspin*rand(1)+1); % случайный выбор номера спина
 % периодические граничные условия
 if Ispin==1  Left=s(Nspin);  else  Left=s(Ispin-1);  end;  if Ispin==Nspin  Right=s(1);  else  Right=s(Ispin+1);  end;  de=2*s(Ispin)*(-h+J*(Left+Right)); % пробное изменение энергии спина
 if de<=Edemon % принятие пробного изменения
 s(Ispin)=-s(Ispin);  Accept=Accept+1;  Edemon=Edemon-de;  Esystem=Esystem+de;  end;  k=k+1;  Es(k)=Esystem;  Ed(k)=Edemon;  SpM(k)=sum(s);  end; end; Accept=Accept/(NTrial*Nspin);

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

>> Nspin=200; % число спинов системы
>> J=1; % константа обменного взаимодействия
>> h=0; % внешнее поле
>> Esi=-10; % задание конечного значения энергии системы
>> NTrial=100; % число испытаний
% вычисление мгновенных значений полной энергии системы (Es),
% энергии демона (Ed), мгновенной намагниченности (SpM), числа
% принятия решений
>> [Es Ed SpM Accept]=Ising(Nspin,J,h,Esi,NTrial);
 % визуализация зависимостей мгновенных значений полной энергии
 % системы и энергии демона от времени
 >> i=1:1200; >> figure(1);plot(i,Es(i), 'k');
>> figure(2);plot(i,Ed(i), 'k'); % вычисление распределения вероятностей мгновенных значений энергии
% демона
>> Nint=50;
>> i=1:Nint;
>> x1=min(Ed);
>> x2=max(Ed);
>> x(i)=x1+(x2-x1)/Nint*i;
>> h=hist(z2,x); % нахождение и визуализация функции, описывающей распределение
% вероятности
>> 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,[10000 0.05]);
>> bar(x,h);colormap white
>> hold on
>> Ni=500;j=1:Ni;
>> X(j)=x1+(x2-x1)/Ni*j;
>> plot(X,F10(beta,X),'k')
>> hold off
>> mean(Es)/Nspin % средняя энергия системы на один спин
ans = -0.0941
>> mean(z3) % средняя намагниченность системы
ans = 1.4352
>> mean(z3)/Nspin % средняя намагниченность на один спин
ans = 0.0072
>> 2/log(1+2./mean(z2)) ans = 9.7904