从简单程序到复杂系统

——用Mathematica探索复杂性


张江

中科院数学与系统科学研究院
复杂系统研究中心


Email : Jakezj@163.com

集智俱乐部：WWW.SWARMAGENTS.COM

“还原论”的局限

现实生活中的例子

小孩喜欢拆卸玩具，却不知道如何组装

科学家碰壁

物理学家的困惑：

经济学家、社会学家的麻烦：人真的是“理性”的吗？

计算机科学家的困惑

“我们人类拆散了整个宇宙，却不知道如何再把它组装起来”
——Barabasi 《连接——网络的新科学》

“复杂”的兴起

《复杂》与圣塔菲研究所（Santa Fe Institute）

Stephen Wolfram也是一名复杂性研究的先行者

复杂性科学的问题：

飞鸟为何聚集成群？

简单、笨拙的蚂蚁为何能够搭建一个等级森严的蚂蚁王国？

无生命的分子如何组成有生命的细胞？

简单的计算机程序如何自发形成网络?

追逐个体利益的自私的人是如何组织构造出秩序井然的市场的？

霍金曾说：“我认为复杂性科学是21世纪的科学”

计算机与复杂系统

计算机游戏的启发

计算机是一个天然的探索复杂系统的工具

计算机可以看作是一个由0、1搭建的复杂系统

20世纪60、70年代，一批科学家在计算机中玩了一些“科学的游戏”，发现了“上帝”的秘密

从一个简单程序开始

Mathematica代码

MyFunc[x_,y_]:={(x^3)-3 x (y^2)+(1/2),3 (x^2) y-(y^3)};
NestJudge[x_,y_]:=Module[{result={x,y},i=0},
For[i=1,i<=10,i++,result=MyFunc[result[[1]],result[[2]]];If[Norm[result]≥10,Break[]]];
If[Abs[result[[1]]]<10||Abs[result[[2]]]<10,Return[1],Return[0]];
];
ArrayPlot[Table[NestJudge[x,y],{x,-5,5,0.05},{y,-5,5,0.05}]]

更多的原始生命

Manipulate[Module[{result},ArrayPlot[Map[(result=NestWhile[BiomorphFunc,#,Norm[#]<R&,1,MaxStep];If[Abs[Re[result]]<R||Abs[Im[result]]<R,1,0])&,Table[x+y I,{x,-Pos[[2]]-radius/2,-Pos[[2]]+radius/2,N[radius/pre]},{y,Pos[[1]]-radius/2,Pos[[1]]+radius/2,N[radius/pre]}],{2}]]],{BiomorphFunc,{Function[z,(z^3)+1/2],Function[z,(z^5)+1/2+I],Function[z,Sin[z]+(z^2)+0.1+1.1 I],Function[z,If[z≠0,(z^z)+(z^3)+3- I,0]],Function[z,If[z≠0,(z^4)+(z^z)+0.1+1.1 I,0]],Function[z,If[z≠0,(z^6)+(z^z)+0.1+I,0]],Function[z,Sin[z]+(E^z)+2+I]},ControlType->PopupMenu},{{Pos,{0,0}},{-2,-2},{2,2}},{{radius,20},1,20,1},{{pre,200},10,1000,10},{{R,10},1,100},{{MaxStep,10},1,100,1},ControlPlacement->Left]

更多的原始生命

Conway的“生命游戏”

1970年，John Conway提出了一个有趣的计算机程序：生命游戏

这是一个方格世界，每个方格仅仅有黑、白两种颜色

“生命游戏”的规则

想象一个外星空间存活着很多方格生命，如果一个方格是黑色的，那么这里就有一个生命，如果是白色的就没有

每个方格周围都有8个邻居

诞生：

如果某时刻，一个空方格 （白色） 周围邻居中有三个活的生命 （黑色方格） ，那么该方格就会诞生一个新生生命 （白色 -> 黑色）

死亡：

在某时刻，对于已经占领某方格的生命体 （黑色） 只有当它的邻居中有2个或者3个是活的 （黑色） ，它才会继续存活，否则不是因为过分拥挤就是因为过分孤独而死亡 （黑色 -> 白色）

实现“生命游戏”

Hiroki Sayama' s Demonstration on demonstration.wolfram.com

Manipulate[
GameRule[lst_]:=If[lst[[2,2]]==0,If[Total[Flatten[lst]]==lp,1,0],If[Total[Flatten[lst]]>4||Total[Flatten[lst]]<3,0,1]];
LifeOneStep:=(Map[GameRule[#]&,RotateRight[RotateRight[#,1]&/@Partition[#,{3,3},{1,1},{1,1}],1],{2}])&;
Refresh[
If[oldlen≠len||oldp≠p,
oldlen=len;oldp=p;updating=False;config=Table[If[RandomReal[]<p,1,0],{len},{len}];t=0
];
If[updating||onestep,
t++;
onestep=False;
config = LifeOneStep[config];
If[config=={},config={{0}}]
];
ArrayPlot[config,PlotLabel->Row[{Style["t",Italic]," = "<>ToString[t]}],ImageSize->{400,400},AspectRatio->Automatic,Mesh->IsMeshed]
,UpdateInterval->If[updating,0,Infinity]
]
,
{{updating,False,"run simulation"},{True,False}},
Button["reset",updating=False;config=Table[If[RandomReal[]<p,1,0],{len},{len}];t=0,ImageSize->Small],
{{bg,False,"infinite space"},{True,False}},
{{len,50,"edge length"},30,200,10,Appearance->"Labeled",ImageSize->Tiny},
{{oldlen,50},ControlType->None},
{{p,0.1,"initial density"},0,1,.01,Appearance->"Labeled",ImageSize->Tiny},
{{oldp,0.1},ControlType->None},
{{t,0},ControlType->None},
{{config,Table[If[RandomReal[]<p,1,0],{50},{50}]},ControlType->None},
{{onestep,False},ControlType->None},{{lp,3,"Updating Rule"},1,9,1,Appearance->"Labeled",ImageSize->Tiny},{{IsMeshed,True},{True,False}},
ControlPlacement->Left,
TrackedSymbols->{updating,onestep,config,len,t,bg,p},AutorunSequencing->{3,5},SynchronousUpdating->True
]

从“生命游戏”到“人工生命”

游戏规则：宇宙的物理法则

生命游戏是一个自己演化的宇宙

Conway曾经自豪的说，如果给我足够大的空间，等待足够长的时间，自然界中的一切事物都能在生命游戏中诞生，包括能够撰写PhD论文的智慧生命

在一个0 、1 数字构成的虚拟世界中，是否存在着数字生命呢？

朗顿的故事

《数字创世纪》 ——李建会、张江

人工生命Boid

靠近、对齐、避免碰撞三条规则涌现出的灵活、生动的复杂行为

Boid被广泛的应用于计算机动画、游戏等领域。

建模Boid

每只Boid都有一个邻居范围

靠近

对齐

避免碰撞

用Mathematica实现Boid

计算机模拟的一般方法

把计算机想象为一个“虚拟宇宙”

内存就是这个“虚拟宇宙”的空间

CPU每“滴答”一次就是“虚拟宇宙”运行一秒

程序代码就是生命体

第一步：系统表示（建立模型）

每一只Boid的位置可以用{x,y}表示，速度用{vx,vy}表示

50只Boid就用50个表表示：

{{{(Subscript[x, 1]),(Subscript[y, 1])}{(Subscript[vx, 1]),(Subscript[vy, 1])}},{{(Subscript[x, 2]),(Subscript[y, 2])}{(Subscript[vx, 2]),(Subscript[vy, 2])}},{{(Subscript[x, 3]),(Subscript[y, 3])}{(Subscript[vx, 3]),(Subscript[vy, 3])}},...}

随机生成50个Boid的坐标和速度：

Dynamic[States]

States=Table[{{RandomReal[{0,100}],RandomReal[{0,100}]},2 {Cos[#],Sin[#]}&@@{RandomReal[{-2 Pi, 2 Pi}]}},{50}]

第二步：将Boid画到屏幕上

Mathematica的画图很简单，只要用Graphics命令，不用考虑坐标转换，例如画一个点

Graphics[Point[{0,0}]]

将50个Boid画到屏幕上：

Graphics[Point[#[[1]]]&/@States]

考虑速度，每个Boid是一个小箭头，写个函数输入{x,y,vx,vy}，画出这个Boid

DrawBoid[state_,Size_,VR_]:=Return[{Arrowheads[0.02],Arrow[{state[[1]],state[[1]]+Size/(VR*20)state[[2]]},0.01]}]

试一下（两只Boid）

Graphics[{DrawBoid[{{1,1},{0.5,0.4}},100,2],DrawBoid[{{0,1},{0.1,0.5}},100,2]}]

把50只画出来

Graphics[DrawBoid[#,100,2]&/@States]

第三步：让Boids动起来

每一次CPU跳动，程序需要做的：

位置 + 速度 = 新的位置，即x = x + vx; y = y + vy

OneStep[lastState_]:={Mod[#[[1]]+#[[2]]+{100,100},100],#[[2]]}&/@lastState

States=OneStep[States];Graphics[DrawBoid[#,100,2]&/@States]

Module[{last=1},Animate[If[last≠t,last=t;States=OneStep[States]];Graphics[DrawBoid[#,100,2]&/@States,PlotRange -> {{0, 100}, {0, 100}}],{t,1,Infinity},AnimationRunning->False]]

第四步：考虑相互作用

{,,}

OneStep[lastState_]:=Module[{Neighbors,Collisions,zero,vNew},
(*得到每个Boid的邻居*)
Neighbors=Table[Join[{lastState[[i]]},Select[lastState,EuclideanDistance[First[#],First[lastState[[i]]]]<10&&EuclideanDistance[First[#],First[lastState[[i]]]]>6&]],{i,1,Length[lastState]}];
(*得到每个Boid即将碰撞上的Boid*)
Collisions=Table[Join[{lastState[[i]]},Select[lastState,EuclideanDistance[First[#],First[lastState[[i]]]]<6&&First[#]≠First[lastState[[i]]]&]],{i,1,Length[lastState]}];
(*0 向量，辅助向量*)
zero=Table[0,{2}];
(*计算三种力更新速度*)
vNew=(0.45 If[Length[#]>1,Normalize[Total[Last/@Rest[#]]/(Length[#]-1)],zero] + 0.1 If[Length[#]>1,Normalize[(Total[First/@Rest[#]]/(Length[#]-1))-First[First[#]]],zero])&/@ Neighbors+0.2 Last/@lastState+(0.25 If[Length[#]>1,Normalize[First[First[#]]-Total[First/@Rest[#]]/(Length[#]-1)],zero]&/@Collisions);
(*Print[vNew];*)
(*将速度的模限定在2以内*)
vNew=2 Normalize/@vNew;
Transpose[{Mod[#[[1]]+#[[2]]+{100,100},100]&/@lastState,vNew}]
]

DynamicModule[{last=1,States=Table[{{RandomReal[{0,100}],RandomReal[{0,100}]},2 {Cos[#],Sin[#]}&@@{RandomReal[{-2 Pi, 2 Pi}]}},{50}]},Animate[If[last≠t,last=t;States=OneStep[States]];Graphics[DrawBoid[#,100,2]&/@States,PlotRange -> {{0, 100}, {0, 100}}],{t,1,Infinity},AnimationRunning->False]]

完整代码

(*初始化*)
DynamicModule[{States,TimeStep=1,lastNums=50,Obstacles},
DrawBoid[state_,Size_,VR_]:=Return[{Arrowheads[0.02],Arrow[{state[[1]],state[[1]]+Size/(VR*20)state[[2]]},0.01]}];
OneStep[lastState_,AdhereD_,CollisionD_,VR_,Size_,CohesionStrength_,AlignmentStrength_,CollisionStrength_,InertiaStrength_,ObstacleStrength_,Obstacles_]:=Module[{Neighbors,Collisions,zero,vNew,rObstacle},
(*得到每个Boid的邻居*)
Neighbors=Table[Join[{lastState[[i]]},Select[lastState,EuclideanDistance[First[#],First[lastState[[i]]]]<AdhereD&&EuclideanDistance[First[#],First[lastState[[i]]]]>CollisionD&]],{i,1,Length[lastState]}];
(*得到每个Boid即将碰撞上的Boid*)
Collisions=Table[Join[{lastState[[i]]},Select[lastState,EuclideanDistance[First[#],First[lastState[[i]]]]<CollisionD&&First[#]≠First[lastState[[i]]]&]],{i,1,Length[lastState]}];
(*0 向量，辅助向量*)
zero=Table[0,{2}];
(*计算三种力更新速度*)
vNew=(AlignmentStrength If[Length[#]>1,Normalize[Total[Last/@Rest[#]]/(Length[#]-1)],zero] + CohesionStrength If[Length[#]>1,Normalize[(Total[First/@Rest[#]]/(Length[#]-1))-First[First[#]]],zero])&/@ Neighbors+InertiaStrength Last/@lastState+(CollisionStrength If[Length[#]>1,Normalize[First[First[#]]-Total[First/@Rest[#]]/(Length[#]-1)],zero]&/@Collisions)+(ObstacleStrength Normalize[Sum[rObstacle=EuclideanDistance[Obstacles[[i]],First[#]];If[rObstacle≠0&&rObstacle<2AdhereD,(First[#]-Obstacles[[i]])/(((EuclideanDistance[Obstacles[[i]],First[#]]))^3),zero],{i,1,Length[Obstacles]}]]&/@lastState);
(*将速度的模限定在VR以内*)
vNew=VR Normalize/@vNew;
Return[Transpose[{Mod[#[[1]]+#[[2]]+{Size,Size},Size]&/@lastState,vNew}]];
];
Manipulate[
If[TimeStep==1,
States=Table[{{RandomReal[{0,Size}],RandomReal[{0,Size}]},VR ({Cos[#],Sin[#]})&@@{RandomReal[{-2 Pi, 2 Pi}]}},{Nums}];
TimeStep++;
];
If[lastNums≠Nums,
If[Nums<lastNums,
States=Drop[States,IntegerPart[lastNums-Nums]];
,
States=Join[States,Table[{{RandomReal[{0,Size}],RandomReal[{0,Size}]},VR ({Cos[#],Sin[#]})&@@{RandomReal[{-2 Pi, 2 Pi}]}},{Nums-lastNums}]];
];
lastNums=Nums;
];
If[!HaveObstacles,Obstacles=Table[{-50,-50},{3}],Obstacles=CObstacles];
If[IsRun&&Mod[TimeStep,Speed]==0,TimeStep++;States=OneStep[States,CohesionD,CollisionD,VR,Size,CohesionStrength,AlignmentStrength,CollisionStrength,InertiaStrength,ObstacleStrength,Obstacles],If[IsRun,TimeStep++;]];Graphics[Join[DrawBoid[#,Size,VR]&/@States,Disk[#,CohesionD/2]&/@Obstacles],PlotRange -> {{0, Size}, {0, Size}}]
,{{IsRun,False},{True,False}},{{Nums,50},1,100,1},{{Size,100},100,1000},{{VR,2},1,500},{{CObstacles,{{0,0},{0,0},{0,0}}},Locator,Visible->HaveObstacles},{{HaveObstacles,False},{True,False}},{{CohesionD,10},4,100},{{CollisionD,6},5,20},{{CohesionStrength,0.1},0,1,0.1},{{AlignmentStrength,0.45},0,1,0.1},{{CollisionStrength,0.2},0,1,0.1},{{InertiaStrength,0.5},0,1,0.1},{{ObstacleStrength,0.7},0,1,0.1},{{Speed,1},1,100,1},ControlPlacement->Left,SynchronousUpdating->True]]

计算机程序能进化吗？

答案是肯定的

遗传算法（Genetic Algorithm）就是模拟自然进化的一种实用算法

将遗传变异、自然选择的规则用计算机实现

《盲目的钟表匠》

由计算机代码生成若干生物体形态

人充当上帝的角色，对这些生物体进行选择

进化可以重演

盲目的钟表匠

盲目的钟表匠源程序

DrawLines[Lines_]:=Graphics[Line[{Lines,{-#[[1]],#[[2]]}&/@Lines}]];
MyTree[PointStart_,order_,dir_,dxy_,outin_]:=Module[{PointNew,NBVAR=8,out,ddir},
ddir=dir;
If[dir<0,ddir=dir+NBVAR,If[dir≥NBVAR,ddir=NBVAR-dir]];
PointNew=PointStart+order*dxy[[ddir+1]];
out=Append[outin,PointNew];
If[order>0,
out=MyTree[PointNew,order-1,ddir-1,dxy,out];
out=MyTree[PointNew,order-1,ddir+1,dxy,out];
];
Return[out];
];
Mutate[Gene_]:=Module[{r1,r2,out,temp},temp=Flatten[Gene];r1=RandomInteger[{1,Length[temp]}];r2=RandomInteger[{-3,3}];out=temp;out[[r1]]=Mod[out[[r1]]+r2,7];Return[Partition[out,2]];]
RandomInitial[]:=Table[{RandomInteger[{-6,6}],RandomInteger[{-6,6}]},{8}]

小试验

Gene={{-2,1},{5,-6},{2,3},{-5,-1},{-1,2},{3,-3},{-1,4},{-4,0}};

DrawLines[MyTree[{0,0},7,8,Gene,{}]]

Gene=Mutate[Gene]

DrawLines[MyTree[{0,0},7,8,Gene,{}]]

交互式小程序

DynamicModule[{Num=12,s},s=Table[RandomInitial[],{Num}];Dynamic[GraphicsGrid[Partition[Table[With[{i=i},Button[DrawLines[MyTree[{0,0},7,0,s[[i]],{}]],s=Table[Mutate[s[[i]]],{Num}],Appearance->None]],{i,Num}],4],Frame->All,ImageSize->{800,800}]]]

共同特征

简单的规则蕴藏着复杂的表现

灵活多变

看待计算机程序的全新视角：模拟的宇宙

复杂性的产生来源于代码之间的相互作用

什么是复杂系统？

涌现— 复杂系统的核心思想

什么是涌现(Emergence)？

简单的个体、微观规律 宏观的复杂行为

整体大于部分之和

该词的含义突然出现的意思，所以涌现也叫突现

涌现的例子

闪烁的霓虹灯

蚂蚁王国

计算机中的涌现

代码的局部相互作用

让规则处于混沌与秩序的边缘