Particle Life粒子生命演化的MATLAB模拟

Particle Life粒子生命演化的MATLAB模拟

  • 0 前言
  • 1 基本原理
    • 1.1 力影响-吸引排斥行为
    • 1.2 距离rmax影响
  • 2 多种粒子相互作用
    • 2.1 双种粒子作用
    • 2.1 多种粒子作用
  • 3 代码

惯例声明:本人没有相关的工程应用经验,只是纯粹对相关算法感兴趣才写此博客。所以如果有错误,欢迎在评论区指正,不胜感激。本文主要关注于算法的实现,对于实际应用等问题本人没有任何经验,所以也不再涉及。

请添加图片描述

0 前言

Particle Life粒子生命演化最早是2017年由数字艺术家Jeffery Ventrella定义的,通过非常简单方法的定义粒子间的作用力,从而产生非常复杂的变化。

最开始Jeffery Ventrella管这种生成方法叫做Clusters,其思想来源于生物学家Lynn Margulus。每个粒子具有不同的颜色,每个颜色代表一种属性。粒子不仅会受到自己颜色粒子的吸引或排斥,也会受到其它颜色粒子的吸引和排斥。

在不同的参数下,粒子间会发生复杂的相互运动,某些参数会呈现出复杂的固定斑图,某些参数会呈现出类似生物之间的集群、逃跑、捕食等各种行为。

章节安排为:第一章主要是讲解原理,第二章演示一些基本的例子,第三章给出了基于MATLAB的具体代码。

本文的参考文献如下:
[1]粒子生命演化:由数量庞大的单体粒子演化出复杂的群体行为逻辑
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 基本原理

首先,假设一群粒子A,它们互相会受到其它粒子的作用力。两个粒子间的力大小是粒子间距离r的函数。

请添加图片描述

当距离r较小,小于rmin时,设置了-1的排斥力,为防止粒子之间重合。当粒子距离在rmin和rmax之间,粒子最大作用力为Fi。当粒子距离超过rmax,设置作用力为0,防止计算量过大。

当然有几个细节点需要注意:

1粒子所受的作用力只遵循上面的力方程,但不一定遵循牛顿第三定理。粒子的速度和加速度通过牛二律F=ma得到。由于防止粒子运动过快,还需要在全场设置粘滞阻尼。所以其实牛顿第一定理也不满足。当然由于这并不是精准的模拟仿真,所以这些小事可以忽略。

2力Fi是可以自行设置的,当Fi<0,粒子间呈现出排斥性,当Fi>0,粒子间呈现出吸引性,一般不超过±2;

3距离rmin通常在rmax的1/4~1/5左右;rmax和画布大小有关,rmax越大,越会有全局的粒子参与,rmax越小,粒子的行为越局部。

1.1 力影响-吸引排斥行为

当F<0时,粒子间呈现出排斥的现象:
请添加图片描述
当F>0时,粒子间呈现出吸引的现象:
请添加图片描述

1.2 距离rmax影响

这里画布大小都定义为1。
当rmax=0.2时,粒子的汇集效果如下:
请添加图片描述
当rmax=0.5时,粒子的汇集效果更全局化:
请添加图片描述

2 多种粒子相互作用

2.1 双种粒子作用

对于两种粒子A和B,力Fi共有4个,分别为A对A之间的力,A对B之间的力,B对A之间的力和B对B之间的力。这4个力可以写为一个矩阵形式:

AB
AF_AAF_AB
BF_BAF_BB

当假设A对A存在吸引,且A还会吸引B。但是B没有反向作用A的力,B与B之间也不会互相作用。这里的矩阵可以写作:
[ 1 0 0.5 0 ] \begin{bmatrix} 1 &0 \\ 0.5&0 \end{bmatrix} [10.500]
此时得到的图形为细胞图案,A粒子在中间互相吸引到一团,周围吸引一圈B粒子。
请添加图片描述
再添加两个规则给粒子B,粒子B之间会弱吸引,但粒子B排斥粒子A。此时由于粒子AB间一个吸引一个排斥,构成了不断向前运动的追逐系统。
[ 1 − 1 0.5 0.5 ] \begin{bmatrix} 1 &-1 \\ 0.5&0.5 \end{bmatrix} [10.510.5]
追逐模型如下:
请添加图片描述
之后多种粒子之间的运动规律,也是由上述各个规则叠加演化而成。
但是由于规则数量等于粒子种类N的平方,比如3种粒子就有9种粒子间规则,4种粒子就有16种粒子间规则。这就导致复杂性暴增,产生了无穷多的变化。

2.1 多种粒子作用

由于规则的复杂性,每一次随机出的结果可能都是独一无二的,且是其它人都未曾见过的。这种随机性和复杂性正是Particle Life的迷人之处。

下面列举一些演示计算结果
三种粒子,细胞图案:
请添加图片描述
三种粒子,岛屿图案:
请添加图片描述

三种粒子,循环捕食图案:
请添加图片描述
5种粒子的交互作用,呈现出一定的结构:
请添加图片描述

3 代码

上面绘图代码见文末。

主要更改粒子数量N,颜色数量NColor即可。建议粒子数量N大概是500倍颜色数量。不易太多,由于MATLAB运行效率较低,所以按照实际电脑配置自行更改。

力的作用距离Rmax在最好是1/c的形式,c是一个整数。

迭代总步数StepMax越大,展示的时间越长。这个如果想长时间欣赏粒子间作用,可以选择一个比较大的数。

图像刷新频率FrameFreq是用来控制多少个时间步显示一次。一般选择2就行,太大会有卡顿的感觉。

clear
clc
close all
%Particle Life粒子生命 MATLAB代码%% 初始设定参数
%初始设定
rng('shuffle');%随机种子
N=1500;%粒子数量
NColor=3;%颜色数量
Ni=rand(NColor,1);Ni=round(Ni*N/sum(Ni));%随机分配每个颜色对应的粒子数量
N=sum(Ni);Rmax=1/5;%力作用的距离
mcp=hsv(40);colormap(mcp(1:32,:));%定义展示颜色
StepMax=1.2e3;%结束迭代时间步
FrameFreq=2;%刷新率,正整数,最小为1,越大图像刷新越慢
%% 其它默认参数
%绘图范围
Xlim=[0,1];
Ylim=[0,1];
%定义每个粒子颜色编号
ColorP=zeros(N,1);
for t=1:NColorColorP(1+sum(Ni(1:t-1)):sum(Ni(1:t)))=t;
end
%粒子的力关系矩阵
FMat=rand(NColor,NColor)*3-1.5;%所有力Fi在-1.5~1.5之间
%粒子坐标速度
XY_P=rand(N,2)*0.8+0.1;%所有粒子点坐标
VXY_P=zeros(N,2);%粒子点速度Rmin=Rmax/5;%粒子间的最小作用距离
MeshMax=1/Rmax;%网格数量
dt=5e-3;%时间精度%构建力函数
t=0;%初始时间
c=Rmax*15.0*sqrt(N);%阻尼,为了防止粒子运动速度太快%% 循环计算每一步迭代
tJ=0;%绘图计数
for kt=1:StepMax%计算点对应的网格XYindx=ceil(XY_P/Rmax);%循环计算每个点所受的力ForceP=zeros(N,2);for kp=1:N %循环每一个点%该点的颜色、坐标和网格Color_k=ColorP(kp,:);XY_k=XY_P(kp,:);XYindx_k=XYindx(kp,:);%计算周围点对该点的力F_k=FMat(Color_k,ColorP)';[Indx_t,XY_P_B,F_B]=Beside9(XYindx_k,XYindx,MeshMax,XY_P,F_k);%周边点索引ForceP_k=F_Func(XY_P_B-XY_k,F_B,Rmin,Rmax);ForceP(kp,:)=ForceP_k;end%增加阻尼项,和v相反ForceP=ForceP-c.*VXY_P;%根据F更新位移x和速度v。dv=at,dx=vt+at^2/2VXY_P_New=VXY_P+ForceP*dt;XY_P=XY_P+0.5*(VXY_P+VXY_P_New)*dt;VXY_P=VXY_P_New;%循环边界条件,如果超出边界,就移到另一端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;%加一时间步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 offpause(0.01)%每一帧图像停留时间tJ=tJ+1;end
end%% 后置函数
function Ft2=F_Func(xy,F,rmin,rmax)
%粒子左右函数
%xy,N行2列的向量,代表别的点距离O点的距离向量
%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;
%计算力向量
dir_xy=xy./r;
dir_xy(isnan(dir_xy))=0;
Ft_Vec=dir_xy.*(Ft*ones(1,2));
%计算合力
Ft2=sum(Ft_Vec,1);
endfunction [BesideIndx1,XY_P_B,F_P]=Beside9(XYindx0,XYindx1,NMesh,XY_P,F_P)
%寻找点0附近区域3×3共9格区域内
%开启循环边界条件%复制出边界点,然后再计算。因为有的点在rmax较大的循环边界条件,会同时向上和下吸引
if XYindx0(1)==1%把最后一列复制一份到前面indx_t=XYindx1(:,1)==NMesh;XYindx1_t=XYindx1(indx_t,:);XYindx1_t(:,1)=0;%赋值为0XYindx1=[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%把第一列复制一份到最后indx_t=XYindx1(:,1)==1;XYindx1_t=XYindx1(indx_t,:);XYindx1_t(:,1)=NMesh+1;%赋值为NMesh+1XYindx1=[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%把最后一行复制一份到前面indx_t=XYindx1(:,2)==NMesh;XYindx1_t=XYindx1(indx_t,:);XYindx1_t(:,2)=0;%赋值为0XYindx1=[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%把第一行复制一份到最后indx_t=XYindx1(:,2)==1;XYindx1_t=XYindx1(indx_t,:);XYindx1_t(:,2)=NMesh+1;%赋值为NMesh+1XYindx1=[XYindx1;XYindx1_t];XY_P=[XY_P;XY_P(indx_t,:)+[0,1]];F_P=[F_P;F_P(indx_t)];
end
%夹在范围之内的点有哪些
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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/117946.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

2023-09-02 LeetCode每日一题(最多可以摧毁的敌人城堡数目)

2023-09-02每日一题 一、题目编号 2511. 最多可以摧毁的敌人城堡数目二、题目链接 点击跳转到题目位置 三、题目描述 给你一个长度为 n &#xff0c;下标从 0 开始的整数数组 forts &#xff0c;表示一些城堡。forts[i] 可以是 -1 &#xff0c;0 或者 1 &#xff0c;其中&…

leetcode 189. 轮转数组

2023.9.3 k的取值范围为0~100000&#xff0c;此时需要考虑到两种情况&#xff0c;当k为0时&#xff0c;此时数组不需要轮转&#xff0c;因此直接return返回&#xff1b;当k大于等于数组nums的大小时&#xff0c;数组将会转为原来的数组&#xff0c;然后再接着轮转&#xff0c;此…

快速上手GIT命令,现学也能登堂入室

系列文章目录 手把手教你安装Git&#xff0c;萌新迈向专业的必备一步 GIT命令只会抄却不理解&#xff1f;看完原理才能事半功倍&#xff01; 快速上手GIT命令&#xff0c;现学也能登堂入室 系列文章目录一、GIT HELP1. 命令文档2. 简要说明 二、配置1. 配置列表2. 增删改查3. …

flutter自定义按钮-文本按钮

目录 前言 需求 实现 前言 最近闲着无聊学习了flutter的一下知识&#xff0c;发现flutter和安卓之间&#xff0c;页面开发的方式还是有较大的差异的&#xff0c;众所周知&#xff0c;android的页面开发都是写在xml文件中的&#xff0c;而flutter直接写在代码里&#xff08;da…

C#搭建WebSocket服务实现通讯

在学习使用websocket之前我们先了解一下websocket&#xff1a; WebSocket是一种在单个TCP连接上进行全双工通信的通信协议。与HTTP协议不同&#xff0c;它允许服务器主动向客户端发送数据&#xff0c;而不需要客户端明确地请求。这使得WebSocket非常适合需要实时或持续通信的应…

分页功能实现

大家好 , 我是苏麟 , 今天聊一聊分页功能 . Page分页构造器是mybatisplus包中的一个分页类 . Page分页 引入依赖 <dependency><groupId>com.baomidou</groupId><artifactId>mybatis-plus-boot-starter</artifactId><version>3.4.1</ver…

MySQL总复习

目录 登录 显示数据库 创建数据库 删除数据库 使用数据库 创建表 添加数据表数据 查询表 添加数据表多条数据 查询表中某数据 增insert 删delete 改update 查select ​ where like ​编辑 范围查找 order by 聚合函数 count max min sum avg g…

正则表达式练习

(function() {//#region 定义正则表达式// const reg /前端/g;// ------------test-------------// const res reg.test("学java,找黑马");// console.log(res)// ------------exec--------------// const res reg.exec("学好前端&#xff0c;找黑马"…

Flutter 状态管理引子

1、为了更好地了解状态管理&#xff0c;先看看什么是状态。 在类似Flutter这样的响应式编程框架中&#xff0c;我们可以认为U相关的开发就是对数据进行封装&#xff0c;将之转换为具体的U1布局或者组件。借用Flutter官网的一张图&#xff0c;可以把我们在第二部分做的所有开发…

高频面试题:如何分别用三种姿势实现三个线程交替打印0到100

最近面试遇到的一道题&#xff0c;需要三个线程交替打印0-100&#xff0c;当时对多线程并不是很熟悉因此没怎么写出来&#xff0c;网上搜了之后得到现 synchronized wait/notifyAll 实现思路&#xff1a;判断当前打印数字和线程数的取余&#xff0c;不等于当前线程则处于等待…

数据结构 day6

1->xmind 2->递归实现程序&#xff1a;输入一个数&#xff0c;输出该数的每一位

取一个整数各偶数位上的数构成一个新的数字

1 题目 我可太难了&#xff0c;这题我的思路有点复杂&#xff0c;遇到的困难很多&#xff0c;总是值传递搞不清楚&#xff0c;地址传递总是写错。 从低位开始取出一个整数s的各奇数位上的数&#xff0c;剩下的偶数位的数依次构成一个新数t。 例如&#xff1a; 输入s&#xff…

软件架构模式+系统架构

架构模式对比 分层模式 一般信息系统中最常见的4层划分如下&#xff1a; Presentation layer 表示层&#xff08;也就是UI层&#xff09;Application layer 应用层&#xff08;也就是服务层&#xff09;Business logic layer 业务逻辑层&#xff08;也就是领域层&#xff09;…

【C++历险记】面向对象|菱形继承及菱形虚拟继承

个人主页&#xff1a;兜里有颗棉花糖&#x1f4aa; 欢迎 点赞&#x1f44d; 收藏✨ 留言✉ 加关注&#x1f493;本文由 兜里有颗棉花糖 原创 收录于专栏【C之路】&#x1f48c; 本专栏旨在记录C的学习路线&#xff0c;望对大家有所帮助&#x1f647;‍ 希望我们一起努力、成长&…

Python 没有 pip 包问题解决

最近需要搞一个干净的Python,从官网上直接下载解压可用的绿色版&#xff0c;发现无法正常使用PiP 一 官网下载Python https://www.python.org/downloads/ 选择 embeddable package,这种是免安装的包&#xff0c;解压后可以直接使用。 二 配置环境变量 添加环境变量&#xff1a…

【Python数据分析】数据分析之numpy基础

实验环境&#xff1a;建立在Python3的基础之上 numpy提供了一种数据类型&#xff0c;提供了数据分析的运算基础&#xff0c;安装方式 pip install numpy导入numpy到python项目 import numpy as np本文以案例的方式展示numpy的基本语法&#xff0c;没有介绍语法的细枝末节&am…

【混合时变参数系统参数估计算法】使用范数总和正则化和期望最大化的混合时变参数系统参数估计算法(Matlab代码实现)

&#x1f4a5;&#x1f4a5;&#x1f49e;&#x1f49e;欢迎来到本博客❤️❤️&#x1f4a5;&#x1f4a5; &#x1f3c6;博主优势&#xff1a;&#x1f31e;&#x1f31e;&#x1f31e;博客内容尽量做到思维缜密&#xff0c;逻辑清晰&#xff0c;为了方便读者。 ⛳️座右铭&a…

MATLAB中circshift函数转化为C语言

背景 有项目算法使用matlab中circshift函数进行运算&#xff0c;这里需要将转化为C语言&#xff0c;从而模拟算法运行&#xff0c;将算法移植到qt。 MATLAB中circshift简单介绍 circshift是循环移位函数。可以使用于数组和矩阵元素的循环移位。 当A是数组 Bcircshift(A,p);如果…

安全学习DAY20_自动化工具项目武器库介绍

信息打点-自动化工具 文章目录 信息打点-自动化工具本节思维导图&概述 各类红蓝队优秀工具项目集合&#xff1a;All-Defense-Tool 自动化-武器库部署F8x 自动化信息搜集-网络空间AsamF 自动化信息搜集-企查信息ENScan 自动化信息搜集-综合架构-ARL&NemoARL灯塔Nemo_Go …

知识图谱实战应用26-基于知识图谱构建《本草纲目》的中药查询与推荐项目应用

大家好,我是微学AI,今天给大家介绍一下知识图谱实战应用26-基于知识图谱构建《本草纲目》的中药查询与推荐项目应用,本文通过Py2neo连接到知识图谱数据库,系统实现了中药的快速查询、关系分析、智能推荐和知识展示等功能。用户可以输入中药的名称或特征进行查询,系统将从知…