元胞自动机基于元胞自动机模拟晶体生长matlab源码

Posted 2021-06-11 Matlab咨询QQ1575304183

1.元胞

元胞又可称为单元。或基元，是元胞自动机的最基本的组成部分。元胞分布在离散的一维、二维或多维欧几里德空间的晶格点上。
每一个元胞都有一个状态，例如下方中间的元胞状态是有一只小虫子，其他元胞的状态是没有小虫子。但如果虫子移动了，那么就是随着时间的推移，状态的变化产生的结果


元胞有很多种，可以是六边形，也可以是三角形等等，我们可以具体问题具体对待

2.元胞空间

元胞所分布在的空间网点集合就是这里的元胞空间。

3.邻居

在一维元胞自动机中，通常以半径，来确定邻居，距离一个元胞，内的所有元胞均被认为是该元胞的邻居。二维元胞自动机的邻居定义较为复杂，但通常有以下几种形式(我们以最常用的规则四方网格划分为例)。下图中黑色元胞为中心元胞，灰色元胞为其邻居，它们的状态一起来计算中心元胞在下一时刻的状态。

也就是说：元胞下一时刻的状态决定于本身状态和它邻居元胞的状态。

4.规则

规则是一个最重要的一点，决定了元胞算法的好坏。

元胞自动机就是根据规则进行局部元胞间的相互作用而引起全局变化。

% DLA
clc;clear;close all;
S=zeros(400,500); % 生成状态矩阵
S(end,:)=1; % 设置状态矩阵中最下面一行元素等于1
A=1;B=1;X=0.8;
rand('state',0); % 设置随机数的状态数
subplot(121);Ii=imshow(1-S,[]); % 显示状态矩阵
T1=title(['times = 1',', total particle=',num2str(sum(S(:)))],...
    'Fontname','times new roman','fontsize',14); % 显示时间与粒子总数
r=rand(1,500);
subplot(122);P1=plot(sum(S,2)/size(S,2),1:size(S,1),'r');% 绘制各行的密度值曲线
set(gca,'Position',[0.57,0.35,0.33,0.36],'YDir','reverse'); % 设置坐标轴属性
xlim([0,max(sum(S,2)/size(S,2))]); % 设置x轴的范围
ylabel('\\ith','fontname','times new roman','fontsize',14); % y轴标注
xlabel('\\it\\rho','fontname','times new roman','fontsize',14); % x轴标注
title('{\\it\\rho} ({\\ith})','fontname','times new roman','fontsize',14); % 加注图题
set(gcf,'DoubleBuffer','on'); % 设置图形窗口的渲染效果
[L1,L2]=size(S); % 返回状态矩阵的行数L1和列数L2
N=500;H=1; % 初始化参数：粒子总数N和时间参数H
h=150; % 设置截顶高度
scale=0.5; % 设置剪切系数
while N<20000;
    R1=2+round([L1-3]*rand); % 随机产生粒子的坐标
    R2=2+round([L2-4]*rand); % 随机产生粒子的坐标
    flag=0; % 控制循环停止的参数
    while R1<L1&R1>1&R2<L2&R2>1&flag==0; % 验证粒子在状态矩阵内部且粒子未被吸附
        he=S(R1,R2-1)+S(R1,R2+1)+S(R1+1,R2); % 计算左、下和右方位的近邻
        if he>0.5; % 判断近邻中有固定粒子
            S(R1,R2)=1; % 运动粒子被吸附
            flag=1; % 标记粒子已经被吸附
        else
            ra=rand; % 粒子进行随机移动的分量
            rb=rand; % 粒子进行随机移动的分量
            R1=R1+(ra>=0.5)-(ra<0.5); % 计算下一时刻粒子的位置坐标
            R2=R2+(rb>=0.5)-(rb<0.5); % 计算下一时刻粒子的位置坐标
        end
    end
    sS=sum(S,2); % 对行所有元素求和
    Se=find(sS);Se=min(Se); % 找出有粒子的最高一行   
    if Se==[size(S,1)-h]; % 判断高度是否达到截顶高度
        Sx=find(S(Se,:)); % 找出最高点粒子的横坐标
        S=cuth(S,h,Se,Sx,scale); % 切去最高点粒子所在的分支
    end
    set(Ii,'CData',1-S); % 显示状态矩阵
    N=sum(S(:)); % 计算粒子总数
    H=H+1; % 累计时间值
    set(P1,'XData',sum(S,2)/size(S,2)); % 更新密度曲线数据
    set(T1,'string',['times = ',num2str(H),', total particle=',num2str(sum(S(:)))]);% 更新时间和粒子总数
    pause(0.02); % 暂停一下，显示动画效果
end

​

完整代码或仿真咨询QQ1575304183