用元胞自动机模型实现晶体的晶相组织动态演化的模拟，规则和算法如下：
1.随机生成一个整数矩阵,矩阵中的每个元素代表一个元胞，元胞的数值相同的代表一个晶粒，不同数值的元胞之间形成自然的晶界。
2.对矩阵中的元素以概率P0进行随机抽取，具体的做法是：生成一个随机数rand，而P0=exp[-Q/RT],其中Q,R,T 是相关的常数，P0也就是一个确定的常数，在每一个时间步内将rand与P0进行比较，如果rand<P0,那么就对矩阵中的元胞进行抽取，没有选中的元胞的数值保持不变，被选中的元胞的数值就采用以下规则3转变：
3.采用的是Moore型邻居规则和周期性的边界条件，
     Moore型邻居规则如下：
              A1            A2           A3
             A4            A5           A6
              A7            A8           A9
周期性的边界条件的意思如下：

A9          A7         A8       A9    A7
A3          A1         A2       A3    A1
A6          A4         A5       A6    A4
A9          A7         A8       A9    A7
A3          A1         A2       A3    A1
周期性的边界条件的意思就是说，只考虑一个N×N矩阵（蓝色所示），但是要生成一个（N+2）×（N+2）的矩阵（黑色所示），这样就使得原来的N×N的矩阵中的每个元素都有8个邻居，也就是左右相邻，上下相邻。
  具体的转变规则如下：
   （1）如果A5的8个邻居元素（上下左右和斜对角）和A5的数值相同，那么A5的数值就保持不变；
    （2） 如果A2，A4, A6,A8中的任意三个元素的数值相同，例如等于a,那么A5的数值转变成a的概率就是P0,否则再判断A1，A3,A7,A9的的数值，如果A1，A3,  A7, A9中有三个元素的数值相等，例如等于b,那么A5的数值就转变成b的概率就是P0；
    （3）如果以上条件都不满足，那么从A5的这8个邻居中任意选取一个元素，A5的数值就转变为这个元素的值的概率也是P0。
然后将此转变规则扩展到整个矩阵中，实现一个时间步CAS的模拟，输出的图形是转变后矩阵中的元素的灰度值。

clc;clear all;
%初始化参数
N=100;         % 矩阵是一个 N×N的
T=1000;           % 参数值
S=500;            % 元胞的可能取值为1-500
Q=300000;           % 参数值
R=831;               % 参数值
p1=exp(-Q/(R*T));   % 计算概率P0值
kk=1000;           % 循环次数即总的时间步数
A=round(1 +499.*rand(N)); %生成随机整数矩阵499=500-1
Ii=imshow(A,[]);         %显示元胞的灰度值
ti=title('CAS=0', 'Fontsize',14,'Fontname','Times New Roman' );%显示标题

   % 循环部分
for k=1:kk;
   
    if    A1==A2==A3==A4==A6==A7==A8==A9==A5&& (rand<P0);
                           A5=A5;
    elseif   A2==A4==A6||A4==A6==A8|| A2==A6==A8&& (rand<P0);
                           A5=A2;
    elseif   A1==A3==A7||A1==A3==A9||A1==A7==A9 && (rand<P0);
                           A5=A1;
    else      A5 (rand<P0)=B(index);  
               
    end
   A=A5;  
%  定义扩展的矩阵Am
      Am=zeros(N+2);           %生成一个全为0的矩阵Am
      Am(2:end-1,2:end-1)=A;  %将A赋值给Am的中间部分
      Am(1,2:end-1)=A(end,:);   %将A最下面一行补到矩阵Am的上边
      Am(end,2:end-1)=A(1,:);    %将A第一行补到矩阵Am的下面
       Am(2:end-1,1)=A(:,end);  %将A最右边一列补到矩阵Am的最左边
      Am(2:end-1,end)=A(:,1);   %将A的第一列补到矩阵Am的最右边
       Am(1,1)=A(end,end);     %将A的右下角补到Am左上角
      Am(end,end)=A(1,1);      %将A的左上角补到Am右下角
      Am(end,1)=A(1,end);       %将A的右上角补到Am左下角
      Am(1,end)=A(end,1);       %将A的左下角补到Am右上角
    rand('state',0);         %生成一个0-1之间的随机数
A1= Am(1:end-2,1:end-2);   %A5的左上近邻     
A2= Am(2:end-1,1:end-2);    %A5的上近邻        
A3 = Am(3:end,1:end-2);      %A5的右上近邻         
A4 = Am(1:end-2,2:end-1);    %A5的左近邻
A6 = Am(3:end,2:end-1);      %A5的右近邻         
A7 = Am(1:end-2,3:end);       %A5的上近邻      
A8 = Am(2:end-1,3:end);        %A5的下近邻     
A9= Am(3:end,3:end);         %A5的右下近邻
A5=Am(end-1,end-1);        %Am的中间部分A5
B={A1;A2;A3;A4;A6;A7;A8;A9};  %将A5的邻居存储在矩阵B中
     index=round(1+7* rand(1));    %生成1-8的参数值便于提取
  set(Ii, 'CData',[]);                %更新灰度值的显示
  set(ti,'string',['CAS=',num2str(k)]);  %显示更新的CAS时间步数
  pause(1);                        %暂停1秒，显示动画效果
end