%-----------------------------------------------------------
function v1=vgb(s,a,n)
v1(1,1)=(a(s(1,1),s(2,1))+a(s(1,1),s(1,2))+a(s(1,1),s(2,2)))/3;
v1(1,n)=(a(s(1,n),s(2,n))+a(s(1,n),s(1,n-1))+a(s(1,n),s(2,n-1)))/3;
v1(n,1)=(a(s(n,1),s(n,2))+a(s(n,1),s(n-1,1))+a(s(n,1),s(n-1,2)))/3;
v1(n,n)=(a(s(n,n),s(n,n-1))+a(s(n,n),s(n-1,n))+a(s(n,n),s(n-1,n-1)))/3;
j=1;
for i=2:n-1
    v1(i,j)=(a(s(i,j),s(i-1,j))+a(s(i,j),s(i+1,j))+a(s(i,j),s(i,j+1))+a(s(i,j),s(i-1,j+1))+a(s(i,j),s(i+1,j+1)))/5;
end
j=n;
for i=2:n-1
v1(i,j)=(a(s(i,j),s(i-1,j))+a(s(i,j),s(i+1,j))+a(s(i,j),s(i,j-1))+a(s(i,j),s(i-1,j-1))+a(s(i,j),s(i+1,j-1)))/5;
end
i=1;
for j=2:n-1
v1(i,j)=(a(s(i,j),s(i+1,j))+a(s(i,j),s(i,j+1))+a(s(i,j),s(i,j-1))+a(s(i,j),s(i+1,j-1))+a(s(i,j),s(i+1,j+1)))/5;
end
i=n;
for j=2:n-1
v1(i,j)=(a(s(i,j),s(i-1,j))+a(s(i,j),s(i,j+1))+a(s(i,j),s(i,j-1))+a(s(i,j),s(i-1,j+1))+a(s(i,j),s(i-1,j-1)))/5;
end
for i=2:n-1
for j=2:n-1
v1(i,j)=(a(s(i,j),s(i-1,j))+a(s(i,j),s(i+1,j))+a(s(i,j),s(i,j+1))+a(s(i,j),s(i,j-1))+a(s(i,j),s(i-1,j-1))+a(s(i,j),s(i-1,j+1))+a(s(i,j),s(i+1,j-1))+a(s(i,j),s(i+1,j+1)))/8;
end
end
%-------------------------------------------------------------
function f=eight1(s,n)
f(1,1)=s(2,1)+s(1,2)+s(2,2);
f(1,n)=s(2,n)+s(1,n-1)+s(2,n-1);
f(n,1)=s(n,2)+s(n-1,1)+s(n-1,2);
f(n,n)=s(n,n-1)+s(n-1,n)+s(n-1,n-1);
j=1;
for i=2:n-1
f(i,j)=s(i-1,j)+s(i+1,j)+s(i,j+1)+s(i-1,j+1)+s(i+1,j+1);
end
j=n;
for i=2:n-1
f(i,j)=s(i-1,j)+s(i+1,j)+s(i,j-1)+s(i-1,j-1)+s(i+1,j-1);
end
i=1;
for j=2:n-1
f(i,j)=s(i+1,j)+s(i,j+1)+s(i,j-1)+s(i+1,j-1)+s(i+1,j+1);
end
i=n;
for j=2:n-1
f(i,j)=s(i-1,j)+s(i,j+1)+s(i,j-1)+s(i-1,j+1)+s(i-1,j-1);
end
for i=2:n-1
for j=2:n-1
    f(i,j)=s(i-1,j)+s(i+1,j)+s(i,j+1)+s(i,j-1)+s(i-1,j-1)+s(i-1,j+1)+s(i+1,j-1)+s(i+1,j+1);
end
end
%-------------------------------------------------------------
function l=eight2(s,n)
l(1,1)=((s(2,1)+s(1,2)+s(2,2))-3)./3;
l(1,n)=((s(2,n)+s(1,n-1)+s(2,n-1))-3)./3;
l(n,1)=((s(n,2)+s(n-1,1)+s(n-1,2))-3)./3;
l(n,n)=((s(n,n-1)+s(n-1,n)+s(n-1,n-1))-3)./3;
j=1;
for i=2:n-1
l(i,j)=((s(i-1,j)+s(i+1,j)+s(i,j+1)+s(i-1,j+1)+s(i+1,j+1))-5)./5;
end
j=n;
for i=2:n-1
l(i,j)=((s(i-1,j)+s(i+1,j)+s(i,j-1)+s(i-1,j-1)+s(i+1,j-1))-5)./5;
end
i=1;
for j=2:n-1
l(i,j)=((s(i+1,j)+s(i,j+1)+s(i,j-1)+s(i+1,j-1)+s(i+1,j+1))-5)./5;
end
i=n;
for j=2:n-1
l(i,j)=((s(i-1,j)+s(i,j+1)+s(i,j-1)+s(i-1,j+1)+s(i-1,j-1))-5)./5;
end
for i=2:n-1
for j=2:n-1
l(i,j)=((s(i-1,j)+s(i+1,j)+s(i,j+1)+s(i,j-1)+s(i-1,j-1)+s(i-1,j+1)+s(i+1,j-1)+s(i+1,j+1))-8)./5;
end
end
%--------------------------------------------------------------
function y1=sum1(s,n)
y1=0;
for i=1:n;
for j=1:n
if s(i,j)==2
y1=y1+1;
end
end
end
%基于CA的产业集群技术发展动力研究
n=20;
time=100;
s=unifrnd(1,2,n,n);
s=round(s);
ran=rand(n);
ran1=rand(n);
ran2=rand(n);
ran3=rand(n);
for t=1:time
a11=0.1;
a12=1;
a21=1.5;
a22=1;
p1=1:-0.01:0;
p2=1:0.02:3;
r=a12*p1(1,t);
m=a21*p2(1,t)
b=[a11,r;m,a22];
w=0.9;
v=zeros(n);
ran=rand(n);
t=1;
vbt{t}=vgb(s,b,n);
st{t}=s;
ft{t}=eight1(s,n);
lt{t}=eight2(s,n);
s;
vbt{t};
ft{t};
lt{t};
for t=1:time;
for i=1:n;
for j=1:n;
s0=st{t};
s00=st{t};
vb00=vgb(s00,b,n);
if s0(i,j)==1
    s0(i,j)=2;
else
s0(i,j)=1;
end
vb0=vgb(s0,b,n);
ft{t}=eight1(st{t},n);
lt{t}=eight2(st{t},n);
jz{t}=w*lt{t};
if vb0(i,j)>=vb00(i,j)&jz{t}(i,j)>0.1&ran1(i,j)<=0.4&ft{t}(i,j)>=3&rand>ran3(i,j);
st{t+1}(i,j)=s0(i,j);
elseif vb0(i,j)>=vb00(i,j)&jz{t}(i,j)>0.1&ran1(i,j)>0.4&ran1(i,j)<=0.7&ft{t}(i,j)>=12&rand>ran3(i,j);
st{t+1}(i,j)=s0(i,j);
elseif vb0(i,j)>=vb00(i,j)&jz{t}(i,j)>0.1&ran1(i,j)>0.7&rand>ran3(i,j);
st{t+1}(i,j)=s00(i,j);
else
st{t+1}(i,j)=s00(i,j);
end
end
end
st{t+1};
vbt{t+1}=vgb(st{t+1},b,n);
vbt{t+1};
ft{t+1}=eight1(st{t+1},n);
ft{t+1};
lt{t+1}=eight2(st{t+1},n);
lt{t+1};
end
end
for t=1:time
y1(t)=sum1(st{t},n);
end
t=1:time;
plot(t,y1)
title('产业集群技术发展动力')
xlabel('模拟次数');ylabel('集群中改进技术的企业总数')