从简单程序到复杂系统
——用Mathematica探索复杂性


北京师范大学系统科学系
张江

“还原论”的局限

现实生活中的例子

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

科学家碰壁

物理学家的困惑：

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

计算机科学家的困惑

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

“复杂”的兴起

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

新的科学问题：

飞鸟为何聚集成群？

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

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

如何创造计算机的智能？

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

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

计算机与复杂系统

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

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

计算机游戏的启发

“给我一个支点，我能撬起整个地球” ——阿基米德

“给我一台计算机，我能创造一个宇宙” ——张江

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

从一个简单程序开始

(x, y) 是平面上的一点，遵循下列的迭代规则

x->(x^3)-3(xy^2)+1/2
y->3 (x^2)y-(y^3)

简单程序的创造

更多的原始生命

用Mathematica实现生物形态

第一步：自定义函数

定义了一个函数叫做MyFunc，该函数的自变量为x,y，输出为新的x'=(x^3)-3 x (y^2)+(1/2)和y'=3 (x^2) y-(y^3)

MyFunc[x_,y_]:={(x^3)-3 x (y^2)+(1/2),3 (x^2) y-(y^3)}

验证一下这个函数

MyFunc[1,0]

MyFunc[3/2,55]

MyFunc[31/8,100]

第二步：迭代

从 {1/2, 1/2} 这个点开始跌代

MyFunc[1/2,1/2]

把{1/4, 1/4} 这个结果再代进MyFunc函数中

MyFunc[1/4,1/4]

N[MyFunc[15/32,1/32]] (*N[]这个函数对括号内容取整*)

用NestList命令从{1/2,1/2}开始自动跌代10次

dlst=NestList[N[MyFunc[#[[1]],#[[2]]]]&,{1/2,1/2},10]

把这些点化到图上（红色），看看它是否超过边界（蓝色的圆）

Show[Plot[Sqrt[100-(x^2)],{x,-10,10}],Plot[-Sqrt[100-(x^2)],{x,-10,10}],ListPlot[dlst,PlotStyle->{Red},Joined->True,Frame->True,PlotMarkers->Automatic],AspectRatio->Automatic,PlotRange->{-10,10}]

第10个点跑出了圆圈，跑出的时候它的x和y都>10，跌代奔跑比较均匀

再从{1,1/2}这个点开始跌代

dlst=NestList[N[MyFunc[#[[1]],#[[2]]]]&,{1,1/2},10]

Show[Plot[Sqrt[100-(x^2)],{x,-10,10}],Plot[-Sqrt[100-(x^2)],{x,-10,10}],ListPlot[dlst,PlotStyle->{Red},Joined->True,Frame->True,PlotMarkers->Automatic],AspectRatio->Automatic,PlotRange->{{-40,10},{-10,10}}]

第4个点跑出了圆圈，跑出的时候|y|<10，只有|x|>10，跌代奔跑不均匀（X方向快）

第三步：迭代+判断

把上面的迭代+判断过程自动化，写成一个自定义函数NestJudge。输入任何一个初始点(x,y)，按照方程MyFunc跌代10次，输出是否奔跑均匀

NestJudge[x_,y_]:=Module[{temp},
temp=NestWhile[N[MyFunc[#1[[1]],#1[[2]]]]&,{x,y},Norm[#1]<10&,1,10];If[Abs[temp[[1]]]<10||Abs[temp[[2]]]<10,Return[1],Return[0]];
]

验证一下：

NestJudge[1/2,1/2]

NestJudge[1,1/2]

第四步：画图

ArrayPlot[Table[NestJudge[x,y],{x,-5,5,0.05},{y,-5,5,0.05}]]

第五步：完整代码

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}]]

Conway的“生命游戏”

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

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

“生命游戏”的规则

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

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

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

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

"Game of Life" in Mathematica

作一个棋盘格世界：

world=Table[RandomInteger[{0,1}],{i,10},{j,10}];

world//MatrixForm

将此世界按照3*3分组：

gworld=RotateRight[RotateRight[#,1]&/@Partition[world,{3,3},{1,1},{1,1}],1];

gworld//MatrixForm

“生命游戏”规则：

GameRule[lst_]:=If[lst[[2,2]]==0,If[Total[Flatten[lst]]==3,1,0],If[Total[Flatten[lst]]>4||Total[Flatten[lst]]<3,0,1]]

nWorld=Map[GameRule[#]&,gworld,{2}];
nWorld//MatrixForm

“生命游戏“迭代：

worlds=NestList[(Map[GameRule[#]&,RotateRight[RotateRight[#,1]&/@Partition[#,{3,3},{1,1},{1,1}],1],{2}])&,Table[RandomInteger[{0,1}],{i,100},{j,100}],10];

ArrayPlot[worlds[[5]],Mesh->True]

探索“生命游戏“：

GameRule[lst_]:=If[lst[[2,2]]==0,If[Total[Flatten[lst]]==3,1,0],If[Total[Flatten[lst]]>4||Total[Flatten[lst]]<3,0,1]];
oneIterate[lst_]:=(Map[GameRule[#]&,RotateRight[RotateRight[#,1]&/@Partition[#,{3,3},{1,1},{1,1}],1],{2}])& [lst];
world=Table[RandomInteger[{0,1}],{i,100},{j,100}];

Manipulate[
 world=oneIterate[world];
ArrayPlot[world,Mesh->True,ImageSize->500],
 {{i, 1, ""}, 1, ∞, 1, ControlType -> Trigger}, 
 SaveDefinitions -> True, AutorunSequencing -> {2, 3, 4}, 
 ControlPlacement -> Left]

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

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

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

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

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

人工生命Boid

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

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

建模Boid

每只Boid都有一个邻居范围

靠近

对齐

避免碰撞

用Mathematica实现Boid

计算机模拟的一般方法

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

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

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

程序代码就是生命体

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

“虚拟空间”：屏幕（长和宽都是100）

每一只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的坐标和速度：

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

第二步：将Boid画到屏幕上

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

Graphics[Point[{0,0}]]

画两个点

Graphics[{Point[{0,0}],Point[{10,10}]}]

将50个Boid画到屏幕上：

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

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

DrawBoid[state_]:=Arrow[{state[[1]],state[[1]]+state[[2]]}]

试一下（两只Boid）

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

把50只画出来

Graphics[DrawBoid/@States]

第三步：让Boids动起来

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

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

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

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

Module[{last=1},Animate[If[last≠t,last=t;States=OneStep[States]];Graphics[DrawBoid[#]&/@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}]
]

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

完整代码

Initialize[]:=Module[{},
vR=2;(*The speed of each boid*)
Size=100;(*The size of the world*)
(*Initialize all positions and velocities of boids*)
rNeighbor=10;(*Neighbor distance*)
rCollision=6;(*Collision distance*)
MySteps=0;(*Simulation steps*)
Init=Table[{{RandomReal[{0,Size}],RandomReal[{0,Size}]},vR {Cos[#],Sin[#]} & @@ RandomReal[{-2*Pi,2*Pi}]},{i,1,50}];
States=Init;(*Simulation states of positions and velocities of each boid in every step*)
];
OneStep[x_List]:=Module[{Neighbors,Collisions,vNew,zero,state},
(*The main loop of the simulation*)
zero=Table[0,{2}];
Neighbors=Table[Join[{x[[i]]},Select[x,EuclideanDistance[First[#],First[x[[i]]]]<rNeighbor&&EuclideanDistance[First[#],First[x[[i]]]]>rCollision&]],{i,1,Length[x]}];
Collisions=Table[Join[{x[[i]]},Select[x,EuclideanDistance[First[#],First[x[[i]]]]<rCollision&&First[#]≠First[x[[i]]]&]],{i,1,Length[x]}];
vNew=(0.3 If[Length[#]>1,Normalize[Total[Last/@Rest[#]]/(Length[#]-1)],zero] + 0.2 If[Length[#]>1,Normalize[(Total[First/@Rest[#]]/(Length[#]-1))-First[First[#]]],zero])&/@ Neighbors+0.2 Last/@x+(0.3 If[Length[#]>1,Normalize[First[First[#]]-Total[First/@Rest[#]]/(Length[#]-1)],zero]&/@Collisions);
vNew=vR Normalize/@vNew;
state=Transpose[{Mod[#+Size,Size]&/@((First[#]+Last[#])&/@x),vNew}];
States=state;
];
DrawIt[x_List]:=Module[{pos,vv},
(*Drawing the boids*)
pos=First[#]&/@ x;
vv=Normalize[Last[#]]&/@x;
{Circle[First[#],1]&/@Transpose[{pos,vv}],Line[{First[#]- Last[#],First[#]+3 Last[#]}]&/@Transpose[{pos,vv}]}
];
ChangeParameters[vv_,num_,sz_]:=Module[{curNum},
(*It will be triggered when parameters changed*)
vR=vv;
Size=sz;
States=Transpose[{Mod[#+Size,Size]&/@First[Transpose[States]],Last[Transpose[States]]}];
curNum=Length[States];
If[curNum>num,States=Drop[States,curNum-num];,
If[curNum<num,States=Join[States,Table[{{RandomReal[{0,Size}],RandomReal[{0,Size}]},vR {Cos[#] , Sin[#]} & @@ RandomReal[{0,2*Pi}]},{i,1,num-curNum}]];
];
];
]
ControlMyAnimate[step_Integer]:=Module[{},
(*Trigger the main simulation loop*)
MySteps=MySteps+1;
OneStep[States];
Return[States];
];
Initialize[]

Manipulate[
 ChangeParameters[vv, num, wsz];
 Graphics[DrawIt[ControlMyAnimate[i]], PlotRange -> {{0, wsz}, {0, wsz}},
   ImageSize -> {ISz, ISz}, 
  Background -> LightYellow], {{ISz, 400, "image size"}, 200, 1000, 
  ImageSize -> Tiny},
 {{vv, 2, "velocity"}, 0.5, 5, ImageSize -> Tiny},
 {{num, 50, "number of Boids"}, 1, 200, 1, ImageSize -> Tiny},
 {{wsz, 100, "world size"}, 50, 1000, 1, ImageSize -> Tiny},
 {{i, 1, ""}, 1, ∞, 1, ControlType -> Trigger}, 
 SaveDefinitions -> True, AutorunSequencing -> {2, 3, 4}, 
 ControlPlacement -> Left]

现实生物的进化

达尔文的进化论

可遗传的变异

优胜劣汰、自然选择

计算机程序能进化吗？

答案是肯定的

遗传算法（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

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

涌现的例子

闪烁的霓虹灯

蚂蚁王国

计算机中的涌现

代码的局部相互作用

Wolfram的《一种新科学》

Stephen Wolfram

15岁开始发表科技论文，22岁获美国MacArthur奖（最年轻的获奖者）

1986年，开创Wolfram研究公司，专门开发Mathematica软件

A New Kind of Science

1980年，他的研究兴趣从物理学转向了复杂性科学

1991年，开始写作巨著《A New Kind of Science》

2002年5月，正式发表厚达1000多页的《A New Kind of Science》

该书的分析与证明没用数学公式，而是用简单程序画出的图形

NKS就是研究各类简单程序的理论与方法

复杂网络

网络无处不在

各类网络有一些相似的性质：小世界、无标度

网络小世界

一个常见的现象：

在一次聚会中，你和一位初次见面的朋友闲聊起来，没多久，你们惊喜地发现，你们都认识某某人，于是感叹： “这个世界真小啊！”

Milgram的“小世界”试验

1961，耶鲁大学的社会心理学家Stanley Milgram做了一个有趣的试验

试验发现，每两个人平均间隔5个人就能发生联系，这个现象被叫做“六度分离”

小世界现象普遍存在于复杂网络中

网络任意两个节点之间的平均距离很小

网络具有聚集效应：你的朋友们也彼此相互熟识

Watts-Strogatz 模型

计算两个指标：平均距离于聚集系数

平均距离：任意两点的最短路径长度的平均值

聚集系数：c(i)=(((∑_(j,k ∈(Subscript[N, i])))(Subscript[a, j,k]))/(Subscript[(C^2), |(Subscript[N, i])|]))

Needs["GraphUtilities`"];
AverageDistance[network_]:=Mean[Mean/@GraphDistanceMatrix[((#/.List->Rule)&/@network)]];ClusteringCoefficient[network_]:=Module[{Neighbors,SNodes,TheMatrix,Clusters,Count=0},
SNodes=Union[First[Transpose[network]],First[Rest[Transpose[network]]]];
Neighbors=Table[DeleteCases[Union[Flatten[Select[network,#[[1]]==SNodes[[i]]||#[[2]]==SNodes[[i]]&]]],SNodes[[i]]],{i,1,Length[SNodes]}];
(*Print[Neighbors];*)
Clusters=(If[Length[#]≤1,0,(Count=0;For[i=1,i≤Length[#]-1,i++,For[j=i+1,j≤Length[#],j++,If[Union[Cases[network,{#[[i]],#[[j]]}],Cases[network,{#[[j]],#[[i]]}]]≠{},Count++;];];];Count
)/(Length[#](Length[#]-1))])&/@Neighbors;
(*Print[Clusters];*)
Return[Mean[Clusters]];
]

从一个规则网络开始

RegNetwork=Flatten[Table[{{i,Mod[i+1,100]},{i,Mod[i+2,100]},{i,Mod[i+3,100]},{i,Mod[i+4,100]}},{i,0,99}],1];
Column[{"Avg Dis:"~~ToString[N[AverageDistance[RegNetwork]]],"Clustering Coef:"~~ToString[N[ClusteringCoefficient[RegNetwork]]],
GraphPlot[#/.{List->Rule}&/@RegNetwork,Method->"CircularEmbedding",ImageSize->500]}]

随机网络

生成一个随机网络

RanNetwork=DeleteCases[Union[Table[{RandomInteger[{0,99}],RandomInteger[{0,99}]},{400}]],{x_,x_}];
Column[{"Avg Dis:"~~ToString[N[AverageDistance[RanNetwork]]],"Clustering Coef:"~~ToString[N[ClusteringCoefficient[RanNetwork]]],
GraphPlot[#/.{List->Rule}&/@RanNetwork,Method->"CircularEmbedding",ImageSize->500]}]

随机改变一些连接

SmallWorldNetworks=Table[Table[
If[RandomReal[{0,1}]<p,{RegNetwork[[i,1]],RandomChoice[DeleteCases[Range[0,99],RegNetwork[[i,1]]]]},RegNetwork[[i]]],{i,Length[RegNetwork]}],{p,0,1,0.01}];

Manipulate[GraphPlot[#/.{List->Rule}&/@SmallWorldNetworks[[IntegerPart[IntegerPart[p*100]]]],Method->"CircularEmbedding",ImageSize->500],{p,0.01,1,0.01}]

p参数对平均距离和聚集系数的影响

result1=Table[dists=DeleteCases[Table[
SmallWorldNetwork=Table[If[RandomReal[{0,1}]<p,{RegNetwork[[i,1]],RandomChoice[DeleteCases[Range[0,99],RegNetwork[[i,1]]]]},RegNetwork[[i]]],{i,Length[RegNetwork]}];AverageDistance[SmallWorldNetwork],{Sample,1,10}],Infinity];{p,N[Mean[dists]]},{p,0,1,0.01}];

result2=Table[Print[p];dists=DeleteCases[Table[
SmallWorldNetwork=Table[If[RandomReal[{0,1}]<p,{RegNetwork[[i,1]],RandomChoice[DeleteCases[Range[0,99],RegNetwork[[i,1]]]]},RegNetwork[[i]]],{i,Length[RegNetwork]}];ClusteringCoefficient[SmallWorldNetwork],{Sample,1,1}],Infinity];{p,N[Mean[dists]]},{p,0,1,0.01}];

TabView[{"Avg Distance"->ListPlot[result1,Frame->True,FrameLabel->{"p","Avg Distance"},Joined->True,PlotMarkers->Automatic,ImageSize->500,AxesOrigin->{0,3}],"Clustering Coefficient"->ListPlot[result2,Frame->True,FrameLabel->{"p","Clustering Coefficient"},Joined->True,PlotMarkers->Automatic,ImageSize->500]}]

集智俱乐部网络社区www.swarmagents.com

集智俱乐部社区网络

从数据库提取网络

file1=#[[1]]&/@Import["/E/mathtalk/from simple to complexity/dbdata1.dat"];
file2=StringSplit[Import["/E/mathtalk/from simple to complexity/dbdata2.dat","Text"]];
file3=#[[1]]&/@Import["/E/mathtalk/from simple to complexity/dbdata3.dat"];
OriDT1={file1,file2,file3};

Poss=Position[OriDT1[[3]],#]&/@Union[DeleteCases[OriDT1[[3]],0.]];
Authors=StringReplace[#,"("~~__->""]&/@OriDT1[[2]];
AuthorGroup=Table[Authors[[#[[1]]]]&/@Poss[[i]],{i,1,Length[Poss]}];
net=Flatten[DeleteCases[ReplaceAll[#,List->Rule]&/@Partition[Union[#],2,1]&/@AuthorGroup,{}]];
netList=net/.Rule->List;
Nodes=Union[Transpose[netList][[1]],Transpose[netList][[2]]];
Print["Number of Nodes:"~~ToString[Length[Nodes]]];
Print["Number of Links:",Length[netList]];

绘制网络

Manipulate[GraphPlot3D[Take[net,i],MultiedgeStyle->False,VertexLabeling->b,Method->MS],{i,1,Length[net],1},{b,{True,False}},{MS,{"SpringEmbedding","SpringElectricalEmbedding","HighDimensionalEmbedding","RandomEmbedding"}}]

平均距离

计算平均距离的函数

Needs["GraphUtilities`"];
AverageDistance[network_]:=Mean[Mean/@#/.List->Rule&/@network];
Print["网络的平均距离",AverageDistance[netList]];
Print["聚集系数",ClusteringCoefficient[netList]];

计算度分布

nLinks=Table[{Nodes[[i]],Select[netList,#[[1]]==Nodes[[i]]||#[[2]]==Nodes[[i]]&]},{i,1,Length[Nodes]}];
degrees=Length[#[[2]]]&/@nLinks;
Tally[degrees];
ListPlot[Tally[degrees],AxesLabel->{"连接数","节点数"}]
fit=Fit[DeleteCases[{Log[#[[1]]],Log[#[[2]]]}&/@Tally[degrees],{___,0}],{1,x},x];
Show[ListLogLogPlot[Tally[degrees],AxesLabel->{"Log[连接数]","Log[节点数]"}],Plot[fit,{x,0,1000},PlotStyle->{Red}]]

按照节点的度的大小顺序列出节点

{#[[1]],Length[#[[2]]]}&/@Sort[nLinks,Length[#1[[2]]]>Length[#2[[2]]]&]//ColumnForm

社团结构

Needs["GraphUtilities`"];
cs=CommunityStructureAssignment[net];
Thread[VertexList[net]->cs];
GraphPlot[net,VertexRenderingFunction->({Hue[cs[[#3]]/Length[Union[cs]]],Point[#]}&),Method->"SpringElectricalEmbedding"]

总结

简单的程序孕育着复杂的行为

复杂性科学主要研究复杂系统的涌现行为

Mathematica是一个强大、便利的探索平台

新的科学研究领域蕴藏着无限的机会