巴克豪森噪声的仿真主要有两种模型:J-A模型和Ising模型,网上参考代码较少,在这里贡献一个自己写的基于伊辛模型的噪声仿真,J-A的有一些问题,感兴趣的可与我联系。
代码主要分三部分:1)磁畴翻转模拟;2)系统哈密顿量计算;3)模型求解。
自旋系统的部分基础参数在代码中,可直接修改,其中sim_time是设置一个蒙特卡洛步数中,磁畴随机翻转的次数,次数不宜太小。通过求解模型,分别画出了初始和结束的系统状态、磁化强度、巴克豪森噪声等图像。
%% 伊辛模型
clear
%% 初始化参数
N=50; %X或Y方向自旋点数
T=3; % 模拟温度
x=round(rand(N,N))*2-1; %初始自旋状态
% x=ones(N,N);
figure(1)
colormap([0 0 0;1 1 1]); %控制两种磁矩的颜色
h=imagesc(x);
title('初始自旋状态图')
MCS=4*10^(3); %Monte Carlo 步数
sim_time=50000;
K_B=0.5; %玻尔兹曼常数 1.3806488*10^(-23)
period=MCS/1000;
H_N=-1;
%归一化外加磁场强度(-1到1之间)
Hn1=-0.998:0.002:1;
Hn2=-Hn1;
for i=1:1:period
if mod(i,2)==1
HN=[H_N Hn1];
H_N=HN;
else
HN=[H_N Hn2];
H_N=HN;
end
end
H_N=H_N';
figure(2)
plot(H_N)
title('归一化外加磁场')
delte=zeros(MCS*sim_time,2);
%% Monte Karlo循环
M=zeros(MCS,1); k=0;
for b=1:1:MCS
for j=1:1:sim_time
Hn=H_N(b);
state_location=randi(N,1,2);
m=state_location(1,1);
n=state_location(1,2);
% Eold=energyf(x,Hn);
x(m,n)=-x(m,n);
% Enew=energyf(x,Hn);
deltaE=energyf_t(x,m,n,Hn);
if deltaE <= 0
x(m,n)=x(m,n);
else
p=rand;
z=exp((-deltaE)/(K_B*T));
if p >= z
x(m,n)=-x(m,n); % 为负显示振荡态
end
end
end
M(b,1)=sum(sum(x))/(N*N);
k=k+1
end
figure(3)
plot(M)
% xlim([2000,MCS])
xlabel('模拟点数');ylabel('幅值')
title('磁化强度')
dm=diff(M);
figure(4)
plot(dm)
xlabel('模拟点数');ylabel('幅值')
title('基于Ising模型的巴克豪森噪声仿真')
figure(5)
colormap([0 0 0;1 1 1]); %控制两种磁矩的颜色
h1=imagesc(x);
title('磁场作用后的自旋状态图')
figure(6)
plot(dm); hold on
plot(H_N(1:end)/4)
xlabel('模拟点数');ylabel('幅值')
title('基于Ising模型的巴克豪森噪声仿真')
%%
figure(7)
plot(H_N(1000:2000),-M(1000:2000));hold on
plot(H_N(2000:3000),-M(2000:3000))
title('磁化强度(归一化形式)')
h1=H_N(1000:2000);h2=H_N(2000:3000);
m1=-M(1000:2000); m2=-M(2000:3000);
y1=smoothdata(m1,'movmean',50); %movmedian'、'gaussian'、'lowess'、'loess'、'rlowess'、'rloess' 或 'sgolay'。
y2=smoothdata(m2,'movmean',50);
figure(8)
plot(h1,y1); hold on; plot(h2,y2)
%%
function deltaE=energyf_t(x,i,j,H_N)
N=50;
J=1;
% 周期性边界条件
if i==1
up=x(N,j);
else
up=x(i-1,j);
end
if i==N
down=x(1,j);
else
down=x(i+1,j);
end
if j==1
left=x(i,N);
else
left=x(i,j-1);
end
if j==N
right=x(i,1);
else
right=x(i,j+1);
end
deltaE=-J*x(i,j).*(up+down+left+right)+H_N*x(i,j);
end
运行结果:
初始系统状态:
4000蒙特卡洛循环后的系统状态:
巴克豪森仿真噪声:
与外加磁场一起显示:
代码实现的效果一般,具体体现在:执行速度慢;巴克豪森仿真噪声峰值之间的区域幅值较大;拓展模型效果差,欢迎大家指正。
另一个注意:通过本代码的系统导出磁化强度画磁滞回线的效果很差,磁滞回线几乎重合。关于伊辛模型的理论知识,可参考北京工业大学王志的博士论文,这篇文章较为详细的写了巴克豪森噪声在伊辛模型下的实现路径,以及模型的拓展。