复杂网络实验9：vicsek模型(matlab)

一、简介

1.定义

为了理解非平衡系统的分簇，传递和相变

vicsek提出了自驱动粒子模型(Self-Propelled Particle,SPP)模型

Vicsek模型是一种用于描述“Active matter”的数学模型

该模型在较高的粒子密度或较低的定位噪声下表现出集体运动

2.数学描述(mathematical description)

3.迭代公式

4.简单理解

刚开始就一群方向随机的粒子

每个周期，每个粒子都检查自己领域半径内粒子的大致方向，得到他们的平均方向

这个方向再加上一个噪声，就是下一时刻的方向

二、实现

1.重要的函数

(1).atan2()

函数名：四象限逆切(four-quadrant inverse tangent)

作用就是给出Y纵坐标和X横坐标，得到角度

再来看看atan()

给定纵坐标，得到角度，但是角度范围在[-pi/2,pi/2]之间，不满足我们要求的360度方向矢量

使用：

拿到周围邻居所有的速度矢量相加，得到邻居的总的矢量，得到Y和X，进而得到总的一个方向趋势

化为角度

(2).quiver()

参考：https://www.cnblogs.com/rgvb178/p/5926168.html

函数名：箭图or矢量图(Quiver or velocity plot)

需要实现这样的效果

#1.出发点为坐标

#2.方向由两个参数 dx 和 dy 决定

    dx为下一时刻横坐标的增量，dy为下一时刻纵坐标的增量

    dx和dy两个矢量相加就是速度矢量

#3.箭头长度不一定为速度矢量的模，可以自己设置

 

实现函数:quiver(x , y , u , v , scale)

x为横坐标的向量，y为纵坐标的向量

u为dx，v为dy，scale为箭头的长度(设置为1箭头长度就为1)

(3).pdist()与squareform()

pdist()方法和squareform()方法一般成对出现

pdist()方法用来计算坐标之间的距离，这里需要计算每个点之间的距离(欧几里得距离)

使用方法：

result = pdist(D,'euclidean') 

D = (x y) x,y均为列向量

D矩阵每一个行向量代表一个点，假设D(m*2)

那么需要计算距离的点为m个，坐标为二维坐标

result是一个向量,这个向量其实是一个矩阵D'(m * m)的下半角矩阵压缩而来的

result其实等于[D'(2,1),D'(2,3),...,D'(2,m),..,D'(m,1),D'(m,m-1)]

result这一向量经过squareform（）函数的反向操作，就能得到这个D'矩阵

D' = squareform(result)

测试：

2.边界重复问题

整个环境是(L*L)的大小，一个粒子越过右边界就要出现在左边界

同样的，越过上边界就要出现在下边界

通过程序的逻辑，我们预期把整个环境做成上下连通，左右连通

最后形成的系统逻辑上是环面(torus)

 

这里，关于边界问题，要考虑两个点

(1).位置迭代的时候矫正位置

这个很好考虑，现在的位置，拿到每个位置的速度矢量，进行位置迭代

拿到新一轮的横坐标和纵坐标，越过边界矫正一下

比如说横坐标范围是[0,L]，有一个粒子这一次迭代得到的横坐标是-1

那么应该加上L，最后经过矫正的横坐标就是 -1 + L

 

(2).计算邻居的时候要考虑到边界的邻居

假如说当前粒子坐标为(0.5,y),邻居搜索半径为1

那么另一头的邻居也要计算到

 

这里只能实现一个粗暴的逻辑，算两遍每个点之间的距离，取相对来说小的那个

首先不要矫正，每个点之间的距离算一遍

然后按照下图的方式，把左边[0,r]范围A内的粒子移动到A'处

同样的把[r,L]范围的粒子移动到左边，再计算一遍每个粒子的距离

中间的粒子坐标就不动

经过这样的方式，在逻辑上满足效果，就是有点浪费性能

3.实现流程

(1).初始化参数

#1.迭代周期

#2.空间大小

#3.噪声大小

#4.粒子数量

#5.作用半径(检测邻居的半径)

#6.速度的模

(2).初始化粒子坐标和角度

(3).进入迭代

#1.画图

#2.越界矫正，得到每个点的距离

#3.得到每个粒子的邻居，计算邻居的平均方向

#4.更新x,y坐标

#5.矫正坐标

#6.更新角度

 

4.全部代码

%%(1).设置初始参数
clear
TMAX = 1000; %迭代1000个周期
L = 7;       %定义长宽为7的一个空间
noise = 0.2;   %噪声大小
N = 200 ;    %粒子数量
r = 1;       %作用半径
vel = 0.03;  %速度的模
%% (2).初始化粒子坐标和角度(1行N列)
x = L * rand(1 , N);
y = L * rand(1 , N);
angle = 2 * pi * rand(1 , N);
%% (3).进入迭代
for step = 1 : TMAX
	% #1.画图
	h = quiver(x , y , cos(angle)  , sin(angle) , 0.3); %箭头长度为0.3
	xlim([0 L]);
	ylim([0 L]);
	axis square
	pause(0.1)
	
    % #2.计算每个点之间的距离
	D = pdist([x' y'] , 'euclidean'); %计算每个点之间的欧氏距离，输入矩阵为N * 2
	%矫正越界坐标，计算相对距离
	tmp_x(x < r) = L + x(x < r);
	tmp_x(x > L-r) = x(x > L-r) - L;
	tmp_x(x >= r & x <= L-r) = x(x >= r & x <= L-r);
	tmp_y(y < r) = L + y(y < r);
	tmp_y(y > L-r) = y(y > L-r) - L;
	tmp_y(y >= r & y <= L-r) = y(y >= r & y <= L-r);
    %计算逻辑矫正的距离
	tmp_D = pdist([tmp_x' tmp_y'] , 'euclidean');
    %两个距离里取到较小的距离
	D = min([D ; tmp_D]); 
	M = squareform(D); %用矩阵形式表示距离
	
    % #3.得到每个粒子的邻居，计算邻居的平均方向
	[l1 , l2] = find(0 < M & M < r); %找到所有距离在r内的距离
	for i = 1:N
		list = l1(l2 == i); %l1记录邻居的下标
		if ~isempty(list)
			new_angle(i) = atan2(mean(sin(angle(list))) , mean(cos(angle(list))));
		else
			new_angle(i) = angle(i);
		end
	end
	% #4.更新x,y坐标
	x = x + vel * cos(angle); 
	y = y + vel * sin(angle);
    
    % #5.世界上下相连，左右相连，矫正坐标
    x(x < 0) = L + x(x < 0);
    x(L < x) = x(L < x) - L;
    y(y < 0) = L + y(y < 0);
    y(L < y) = y(L < y) - L;
    % #6.更新角度
	angle = new_angle + noise * (rand(1,N)); 
    

    
    %窗口关闭，程序退出
    if findobj==0  
     break;
    end
end