巴克豪森噪声的Matlab仿真

打工仔日记

已于 2022-05-30 20:40:33 修改

阅读量701

点赞数 2

文章标签： matlab

于 2022-05-30 10:37:42 首次发布

本文链接：https://blog.csdn.net/qq_46022231/article/details/125042131

版权

巴克豪森噪声的仿真主要有两种模型：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

运行结果：

初始系统状态：