慣例聲明：本人沒(méi)有相關(guān)的工程應(yīng)用經(jīng)驗(yàn)，只是純粹對(duì)相關(guān)算法感興趣才寫此博客。所以如果有錯(cuò)誤，歡迎在評(píng)論區(qū)指正，不勝感激。本文主要關(guān)注于算法的實(shí)現(xiàn)，對(duì)于實(shí)際應(yīng)用等問(wèn)題本人沒(méi)有任何經(jīng)驗(yàn)，所以也不再涉及。

0 前言

Particle Life粒子生命演化最早是2017年由數(shù)字藝術(shù)家Jeffery Ventrella定義的，通過(guò)非常簡(jiǎn)單方法的定義粒子間的作用力，從而產(chǎn)生非常復(fù)雜的變化。

最開(kāi)始Jeffery Ventrella管這種生成方法叫做Clusters，其思想來(lái)源于生物學(xué)家Lynn Margulus。每個(gè)粒子具有不同的顏色，每個(gè)顏色代表一種屬性。粒子不僅會(huì)受到自己顏色粒子的吸引或排斥，也會(huì)受到其它顏色粒子的吸引和排斥。

在不同的參數(shù)下，粒子間會(huì)發(fā)生復(fù)雜的相互運(yùn)動(dòng)，某些參數(shù)會(huì)呈現(xiàn)出復(fù)雜的固定斑圖，某些參數(shù)會(huì)呈現(xiàn)出類似生物之間的集群、逃跑、捕食等各種行為。

章節(jié)安排為：第一章主要是講解原理，第二章演示一些基本的例子，第三章給出了基于MATLAB的具體代碼。

本文的參考文獻(xiàn)如下：
[1]粒子生命演化：由數(shù)量龐大的單體粒子演化出復(fù)雜的群體行為邏輯
https://www.bilibili.com/video/BV1Dh4y1t7hn/
https://www.youtube.com/watch?v=p4YirERTVF0
[2]https://particle-life.com
[3]blender3.6模擬-粒子生命-Particle Life
https://www.bilibili.com/video/BV1Ns4y1B7Fu/

1 基本原理

首先，假設(shè)一群粒子A，它們互相會(huì)受到其它粒子的作用力。兩個(gè)粒子間的力大小是粒子間距離r的函數(shù)。

當(dāng)距離r較小，小于rmin時(shí)，設(shè)置了-1的排斥力，為防止粒子之間重合。當(dāng)粒子距離在rmin和rmax之間，粒子最大作用力為Fi。當(dāng)粒子距離超過(guò)rmax，設(shè)置作用力為0，防止計(jì)算量過(guò)大。

當(dāng)然有幾個(gè)細(xì)節(jié)點(diǎn)需要注意：

1粒子所受的作用力只遵循上面的力方程，但不一定遵循牛頓第三定理。粒子的速度和加速度通過(guò)牛二律F=ma得到。由于防止粒子運(yùn)動(dòng)過(guò)快，還需要在全場(chǎng)設(shè)置粘滯阻尼。所以其實(shí)牛頓第一定理也不滿足。當(dāng)然由于這并不是精準(zhǔn)的模擬仿真，所以這些小事可以忽略。

2力Fi是可以自行設(shè)置的，當(dāng)Fi<0，粒子間呈現(xiàn)出排斥性，當(dāng)Fi>0，粒子間呈現(xiàn)出吸引性，一般不超過(guò)±2；

3距離rmin通常在rmax的1/4~1/5左右；rmax和畫布大小有關(guān)，rmax越大，越會(huì)有全局的粒子參與，rmax越小，粒子的行為越局部。

1.1 力影響-吸引排斥行為

當(dāng)F<0時(shí)，粒子間呈現(xiàn)出排斥的現(xiàn)象：

當(dāng)F>0時(shí)，粒子間呈現(xiàn)出吸引的現(xiàn)象：

1.2 距離rmax影響

這里畫布大小都定義為1。
當(dāng)rmax=0.2時(shí)，粒子的匯集效果如下：

當(dāng)rmax=0.5時(shí)，粒子的匯集效果更全局化：

2 多種粒子相互作用

2.1 雙種粒子作用

對(duì)于兩種粒子A和B，力Fi共有4個(gè)，分別為A對(duì)A之間的力，A對(duì)B之間的力，B對(duì)A之間的力和B對(duì)B之間的力。這4個(gè)力可以寫為一個(gè)矩陣形式：

當(dāng)假設(shè)A對(duì)A存在吸引，且A還會(huì)吸引B。但是B沒(méi)有反向作用A的力，B與B之間也不會(huì)互相作用。這里的矩陣可以寫作：
$$\left[\begin{array}{cc}1& 0\\ 0.5& 0\end{array}\right]$$
此時(shí)得到的圖形為細(xì)胞圖案，A粒子在中間互相吸引到一團(tuán)，周圍吸引一圈B粒子。

再添加兩個(gè)規(guī)則給粒子B，粒子B之間會(huì)弱吸引，但粒子B排斥粒子A。此時(shí)由于粒子AB間一個(gè)吸引一個(gè)排斥，構(gòu)成了不斷向前運(yùn)動(dòng)的追逐系統(tǒng)。
$$\left[\begin{array}{cc}1& ?1\\ 0.5& 0.5\end{array}\right]$$
追逐模型如下：

之后多種粒子之間的運(yùn)動(dòng)規(guī)律，也是由上述各個(gè)規(guī)則疊加演化而成。
但是由于規(guī)則數(shù)量等于粒子種類N的平方，比如3種粒子就有9種粒子間規(guī)則，4種粒子就有16種粒子間規(guī)則。這就導(dǎo)致復(fù)雜性暴增，產(chǎn)生了無(wú)窮多的變化。

2.1 多種粒子作用

由于規(guī)則的復(fù)雜性，每一次隨機(jī)出的結(jié)果可能都是獨(dú)一無(wú)二的，且是其它人都未曾見(jiàn)過(guò)的。這種隨機(jī)性和復(fù)雜性正是Particle Life的迷人之處。

下面列舉一些演示計(jì)算結(jié)果
三種粒子，細(xì)胞圖案：

三種粒子，島嶼圖案：

三種粒子，循環(huán)捕食圖案：

5種粒子的交互作用，呈現(xiàn)出一定的結(jié)構(gòu)：

3 代碼

上面繪圖代碼見(jiàn)文末。

主要更改粒子數(shù)量N，顏色數(shù)量NColor即可。建議粒子數(shù)量N大概是500倍顏色數(shù)量。不易太多，由于MATLAB運(yùn)行效率較低，所以按照實(shí)際電腦配置自行更改。

力的作用距離Rmax在最好是1/c的形式，c是一個(gè)整數(shù)。

迭代總步數(shù)StepMax越大，展示的時(shí)間越長(zhǎng)。這個(gè)如果想長(zhǎng)時(shí)間欣賞粒子間作用，可以選擇一個(gè)比較大的數(shù)。

圖像刷新頻率FrameFreq是用來(lái)控制多少個(gè)時(shí)間步顯示一次。一般選擇2就行，太大會(huì)有卡頓的感覺(jué)。文章來(lái)源地址http://www.zghlxwxcb.cn/news/detail-682823.html

clear
clc
close all
%Particle Life粒子生命 MATLAB代碼

%% 初始設(shè)定參數(shù)
%初始設(shè)定
rng('shuffle');%隨機(jī)種子
N=1500;%粒子數(shù)量
NColor=3;%顏色數(shù)量
Ni=rand(NColor,1);Ni=round(Ni*N/sum(Ni));%隨機(jī)分配每個(gè)顏色對(duì)應(yīng)的粒子數(shù)量
N=sum(Ni);

Rmax=1/5;%力作用的距離
mcp=hsv(40);colormap(mcp(1:32,:));%定義展示顏色
StepMax=1.2e3;%結(jié)束迭代時(shí)間步
FrameFreq=2;%刷新率，正整數(shù)，最小為1，越大圖像刷新越慢
%% 其它默認(rèn)參數(shù)
%繪圖范圍
Xlim=[0,1];
Ylim=[0,1];
%定義每個(gè)粒子顏色編號(hào)
ColorP=zeros(N,1);
for t=1:NColor
    ColorP(1+sum(Ni(1:t-1)):sum(Ni(1:t)))=t;
end
%粒子的力關(guān)系矩陣
FMat=rand(NColor,NColor)*3-1.5;%所有力Fi在-1.5~1.5之間
%粒子坐標(biāo)速度
XY_P=rand(N,2)*0.8+0.1;%所有粒子點(diǎn)坐標(biāo)
VXY_P=zeros(N,2);%粒子點(diǎn)速度


Rmin=Rmax/5;%粒子間的最小作用距離
MeshMax=1/Rmax;%網(wǎng)格數(shù)量
dt=5e-3;%時(shí)間精度

%構(gòu)建力函數(shù)
t=0;%初始時(shí)間
c=Rmax*15.0*sqrt(N);%阻尼，為了防止粒子運(yùn)動(dòng)速度太快

%% 循環(huán)計(jì)算每一步迭代
tJ=0;%繪圖計(jì)數(shù)
for kt=1:StepMax
    %計(jì)算點(diǎn)對(duì)應(yīng)的網(wǎng)格
    XYindx=ceil(XY_P/Rmax);
    %循環(huán)計(jì)算每個(gè)點(diǎn)所受的力
    ForceP=zeros(N,2);
    for kp=1:N %循環(huán)每一個(gè)點(diǎn)
        %該點(diǎn)的顏色、坐標(biāo)和網(wǎng)格
        Color_k=ColorP(kp,:);
        XY_k=XY_P(kp,:);
        XYindx_k=XYindx(kp,:);
        %計(jì)算周圍點(diǎn)對(duì)該點(diǎn)的力
        F_k=FMat(Color_k,ColorP)';
        
        [Indx_t,XY_P_B,F_B]=Beside9(XYindx_k,XYindx,MeshMax,XY_P,F_k);%周邊點(diǎn)索引

        ForceP_k=F_Func(XY_P_B-XY_k,F_B,Rmin,Rmax);
        ForceP(kp,:)=ForceP_k;
    end
    %增加阻尼項(xiàng)，和v相反
    ForceP=ForceP-c.*VXY_P;

    %根據(jù)F更新位移x和速度v。dv=at,dx=vt+at^2/2
    VXY_P_New=VXY_P+ForceP*dt;
    XY_P=XY_P+0.5*(VXY_P+VXY_P_New)*dt;
    VXY_P=VXY_P_New;

    %循環(huán)邊界條件，如果超出邊界，就移到另一端
    XY_P(XY_P>1)=XY_P(XY_P>1)-1;
    XY_P(XY_P<0)=XY_P(XY_P<0)+1;

    t=t+dt;%加一時(shí)間步
    if ~mod(kt,FrameFreq)
        f=figure(1);
        f.Color=[1,1,1];
        cla;
        scatter(XY_P(:,1),XY_P(:,2),6,ColorP,"filled");
        xlim([0,1]);ylim([0,1]);
        %set(gca,'XTick',[],'YTick',[])
        axis off
        pause(0.01)%每一幀圖像停留時(shí)間
        tJ=tJ+1;
    end
end

%% 后置函數(shù)
function Ft2=F_Func(xy,F,rmin,rmax)
%粒子左右函數(shù)
%xy，N行2列的向量，代表別的點(diǎn)距離O點(diǎn)的距離向量
%F，N行1列的向量，代表吸引力F大小
rmid=0.5*(rmax+rmin);
dmid=0.5*(rmax-rmin);
r=sqrt(xy(:,1).^2+xy(:,2).^2);%距離
%r(r==0)=rmax;
Ft=zeros(size(r));
%第一段
indx1=(r<rmin);
Ft(indx1)=r(indx1)/rmin-1;
%第二段
indx_last=~indx1;
indx2=indx_last&(r<rmid);
Ft(indx2)=F(indx2).*(r(indx2)-rmin)/dmid;
%第三段
indx3=(r>=rmid)&(r<rmax);
Ft(indx3)=-F(indx3).*(r(indx3)-rmax)/dmid;
%計(jì)算力向量
dir_xy=xy./r;
dir_xy(isnan(dir_xy))=0;
Ft_Vec=dir_xy.*(Ft*ones(1,2));
%計(jì)算合力
Ft2=sum(Ft_Vec,1);
end

function [BesideIndx1,XY_P_B,F_P]=Beside9(XYindx0,XYindx1,NMesh,XY_P,F_P)
%尋找點(diǎn)0附近區(qū)域3×3共9格區(qū)域內(nèi)
%開(kāi)啟循環(huán)邊界條件

%復(fù)制出邊界點(diǎn)，然后再計(jì)算。因?yàn)橛械狞c(diǎn)在rmax較大的循環(huán)邊界條件，會(huì)同時(shí)向上和下吸引
if XYindx0(1)==1
    %把最后一列復(fù)制一份到前面
    indx_t=XYindx1(:,1)==NMesh;
    XYindx1_t=XYindx1(indx_t,:);
    XYindx1_t(:,1)=0;%賦值為0
    XYindx1=[XYindx1;XYindx1_t];
    XY_P=[XY_P;XY_P(indx_t,:)+[-1,0]];
    F_P=[F_P;F_P(indx_t)];
end
if XYindx0(1)==NMesh
    %把第一列復(fù)制一份到最后
    indx_t=XYindx1(:,1)==1;
    XYindx1_t=XYindx1(indx_t,:);
    XYindx1_t(:,1)=NMesh+1;%賦值為NMesh+1
    XYindx1=[XYindx1;XYindx1_t];
    XY_P=[XY_P;XY_P(indx_t,:)+[1,0]];
    F_P=[F_P;F_P(indx_t)];
end
if XYindx0(2)==1
    %把最后一行復(fù)制一份到前面
    indx_t=XYindx1(:,2)==NMesh;
    XYindx1_t=XYindx1(indx_t,:);
    XYindx1_t(:,2)=0;%賦值為0
    XYindx1=[XYindx1;XYindx1_t];
    XY_P=[XY_P;XY_P(indx_t,:)+[0,-1]];
    F_P=[F_P;F_P(indx_t)];
end
if XYindx0(2)==NMesh
    %把第一行復(fù)制一份到最后
    indx_t=XYindx1(:,2)==1;
    XYindx1_t=XYindx1(indx_t,:);
    XYindx1_t(:,2)=NMesh+1;%賦值為NMesh+1
    XYindx1=[XYindx1;XYindx1_t];
    XY_P=[XY_P;XY_P(indx_t,:)+[0,1]];
    F_P=[F_P;F_P(indx_t)];
end
%夾在范圍之內(nèi)的點(diǎn)有哪些
BesideIndx_X=(XYindx0(1)-1<=XYindx1(:,1))&(XYindx1(:,1)<=XYindx0(1)+1);
BesideIndx_Y=(XYindx0(2)-1<=XYindx1(:,2))&(XYindx1(:,2)<=XYindx0(2)+1);
BesideIndx1=BesideIndx_X & BesideIndx_Y;
XY_P_B=XY_P(BesideIndx1,:);
F_P=F_P(BesideIndx1);
end

到了這里，關(guān)于Particle Life粒子生命演化的MATLAB模擬的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！