【模拟】萤火虫的协同

上万个萤火虫会协同发光，但是其中没有一个领导指挥

那么，其原理何在呢？这里用30个萤火虫进行一下简单模拟

规则：

萤火虫的原始位置、原始状态完全随机
萤火虫的周期也是随机的
某个萤火虫被别的萤火虫影响，如果整体更快（相位大），则被“拖快”；如果整体更慢（相位小），则被“拖慢”。其影响程度受耦合强度 $K$ 影响

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# =====================================
# 1. 参数设置
# =====================================
N = 30  # 萤火虫数量
theta = np.random.rand(N) * 2 * np.pi  # 初始相位，随机在 [0, 2π]
omega = np.random.normal(1.0, 0.05, N)  # 自然频率，有微弱个体差异
K = 1.0  # 耦合强度（影响同步速度）
dt = 0.1  # 时间步长

# 萤火虫的位置
loc = np.random.rand(N, 2) * 10

# =====================================
# 2. 可视化初始化
# =====================================
fig, ax = plt.subplots()
scat = ax.scatter(loc[:, 0], loc[:, 1], s=80, c='black')  # 初始点颜色
ax.set_xlim(0, 10)
ax.set_ylim(0, 10)
ax.set_title('Kuramoto-firefly')
frame_text = ax.text(0.02, 9.5, '', fontsize=12, color='blue')


# =====================================
# 3. 相位 → 颜色映射函数
# =====================================
def phase_to_color(theta_vals):
    norm_vals = (1 + np.cos(theta_vals)) / 2  # 将 cos(θ) 映射到 [0,1]
    return plt.cm.YlOrBr(norm_vals)  # 使用 YlOrBr 配色显示发光强度


# =====================================
# 4. 动画更新函数
# =====================================
def update(frame):
    global theta

    # 计算耦合项（平均相位影响）
    theta_diff = theta[:, np.newaxis] - theta[np.newaxis, :]
    coupling = np.sin(-theta_diff).mean(axis=1)  # 平均耦合项

    # 更新每个萤火虫的相位
    theta += dt * (omega + K * coupling)
    theta = np.mod(theta, 2 * np.pi)  # 保持在 [0, 2π] 范围内

    # 更新图像颜色和文本
    scat.set_color(phase_to_color(theta))
    frame_text.set_text(f'iter: {frame}')
    return scat, frame_text


# =====================================
# 5. 运行动画
# =====================================
ani = FuncAnimation(fig, update, frames=300, interval=100, blit=True)
plt.show()
ani.save("firefly.gif", writer='pillow', fps=10)

数学解释：Kuramoto模型

它是耦合振荡器同步行为的一个经典模型，可以用来解释萤火虫同步发光、青蛙同步叫声、心脏起搏细胞同步跳动等现象。

模型简述：

\[\frac{d\theta_i}{dt} = \omega_i + \frac{K}{N} \sum_{j=1}^N \sin(\theta_j - \theta_i)\]

其中：

$\theta_i$：第 $i$ 个萤火虫的相位（一个周期代表一次发光）
$\omega_i$：其自然频率（本体节奏）
$K$：耦合强度（受他人影响的程度）
$N$：萤火虫总数

关键机制在于“非线性耦合”和“正反馈”：

自发节律： 每只萤火虫都有自己独立的节奏（(\omega_i) 不同）。
局部相互作用： 每只虫只受邻近虫的影响。
正反馈耦合： 如果你发现邻居发光，你更可能调整节奏靠近他们。
非线性同频锁相： 经过若干轮调整，所有虫的节奏渐渐同步。

其它

这篇文章是2015年写的，当时是用 Matlab 写（如下），后来用 Python 重写，并添加了数学描述（如上）

%萤火虫

loc=rand(30,2)*10;

c3=rand(30,1)*10;%能量态

light=ones(30,1);%逻辑态

for i=1:10000   
 scatter(loc(:,1),loc(:,2),light*20+10*ones(30,1))
    c3=c3-light*4+2+sum(light)/20;    %受别的萤火虫影响程度   

for j=1:30%点亮的判断   

if c3(j)<=0      
  light(j)=0;  
  elseif c3(j)>=0     
   light(j)=1;  
   end
end
i   
pause(0.1)
end