数学之美:基于python/matlab实现的赛博数字生命

程序员架构进阶

发布于 2025-11-21 17:45:51

一赛博脊虫与数字生命

最近看到一些视频号发了很多所谓"数字生命"的视频，例如如下的"赛博脊虫"就是其中的一个典型示例。这些数字生命，实际上是数学公式的可视化实现。本篇就以实际代码进行说明，作为探索数字生命、体验数学之美的钥匙。

二实现工具

在本篇之前，我曾经简单探索过Processing，以及manim，但发现对数学公式的可视化实现效果并不是很好。对这类数学公式类的实现，还是matlab或python+matplotlib效果要更好些，无论是在学习成本还是展现效果上。所以本篇基于matlab和python来给出代码示例。

一些基础的工具安装这里不再做详细描述，你可以从matlab官网(https://ww2.mathworks.cn/products/matlab.html)下载matlab并安装，或本地安装python3，并用pip安装numpy和matplotlib的依赖。以下每个图例都会给出两种代码实现。

三赛博脊虫的代码实现与公式分析

话不多说，先放代码。

3.1 matlab实现

% 赛博脊虫
figure('Position',[300,50,900,900], 'Color','k');
axes(gcf, 'NextPlot','add', 'Position',[0,0,1,1], 'Color','k');
axis([0, 400, 0, 400])
SHdl = scatter([], [], 2, 'filled','o','w', 'MarkerEdgeColor','none', 'MarkerFaceAlpha',.4);
t = 0;
i = 0:1e4;
x = i;
y = i./235;
e = y./8 - 13;
while true
    t = t + pi/720;
    k = (4 + sin(y.*2 - t).*3).*cos(x./29);
    d = vecnorm([k; e]);
    q = 3.*sin(k.*2) + 0.3./k + sin(y./25).*k.*(9 + 4.*sin(e.*9 - d.*3 + t.*2));
    SHdl.XData = q + 30.*cos(d - t) + 200; 
    SHdl.YData = 620 - q.*sin(d - t) - d.*39;
    drawnow
end

如下图所示，在matlab中新建文件，把上述代码贴进去，并点击上方的运行按钮，我们就能在弹窗中看到这条赛博脊虫运动的状态了。

3.2 python

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

# 设置图形和坐标轴
fig, ax = plt.subplots(figsize=(9, 9), facecolor='black')
ax.set_facecolor('black')
ax.set_xlim(0, 400)
ax.set_ylim(0, 400)
ax.set_aspect('equal')
ax.axis('off')  # 关闭坐标轴刻度和标签

# 初始化散点对象
scat = ax.scatter([], [], s=2, c='white', alpha=0.4, marker='o', edgecolors='none')

# 预计算一些常量和数据
i = np.arange(10000)  # 使用 10000 个点
x = i.astype(float)
y = i / 235.0
e = y / 8.0 - 13.0

# 动画更新函数
def update(frame):
    t = frame * np.pi / 720  # 计算当前时间 t

    # 计算 k 和 d
    k = (4 + np.sin(y * 2 - t) * 3) * np.cos(x / 29)
    d = np.sqrt(k**2 + e**2)  # 使用 np.linalg.norm([k, e], axis=0) 也可以，但 sqrt 更直接

    # 计算 q
    q = 3 * np.sin(k * 2) + 0.3 / k + np.sin(y / 25) * k * (9 + 4 * np.sin(e * 9 - d * 3 + t * 2))

    # 计算新的 X 和 Y 数据
    new_x = q + 30 * np.cos(d - t) + 200
    new_y = 620 - q * np.sin(d - t) - d * 39

    # 更新散点的数据
    scat.set_offsets(np.column_stack((new_x, new_y)))

    return scat,

# 创建动画
# interval: 每帧之间的间隔（毫秒），值越小动画越快
# blit: 优化绘制，只重绘变化的部分（可能在某些后端不工作）
ani = FuncAnimation(fig, update, interval=20, blit=False, cache_frame_data=False) # blit=False 更兼容

# 显示图形窗口
plt.tight_layout()
plt.show()

# 注意：在某些环境中，plt.show() 会阻塞程序，动画会在此处运行
# 在交互式环境（如 Jupyter Notebook）中，可能需要特殊的后端支持

保存到文件中，例如:cyber_worm.py。

在命令行/终端执行python cyber_worm.py即可。下图左侧是python的实现，右图是matlab的效果。

3.3 公式说明

3.3.1 核心变形函数

包括以下三个：

k = (4 + sin(2y - t) × 3) × cos(x/29) 
d = √(k² + e²)       
q = 3×sin(2k) + 0.3/k + sin(y/25)×k×(9 + 4×sin(9e - 3d + 2t))

详细说明如下：

1、k = (4 + sin(2y - t) × 3) × cos(x/29)

这是一个振幅调制的余弦函数，

振幅部分：4 + 3×sin(2y - t)，随时间t和y坐标变化，

波形部分：cos(x/29)，沿x轴的周期性变化

2、距离参数 d = √(k² + e²)，向量[k,e]的欧几里得范数

3、变形参数 q = 3×sin(2k) + 0.3/k + sin(y/25)×k×(9 + 4×sin(9e - 3d + 2t))

这个复杂函数包含三个部分：

第一部分：3×sin(2k) - 基于k的正弦波动

第二部分：0.3/k - 双曲线分量（当k接近0时产生大值）

第三部分：sin(y/25)×k×(9 + 4×sin(9e - 3d + 2t)) - 复杂的调制波动

3.3.2 最终坐标变换

X坐标变换：X = q + 30×cos(d - t) + 200

基础位置：q

圆形运动：30×cos(d - t)，产生围绕中心的旋转

中心偏移：+200，将图形移到画布中央

Y坐标变换：Y = 620 - q×sin(d - t) - 39×d

基础位置：620（从顶部开始）

圆形运动：-q×sin(d - t)，与X坐标的cos分量配合形成圆形轨迹

径向偏移：-39×d，基于距离d的线性偏移

3.3.3 整体效果分析

本质上是一个参数化的粒子系统，每个粒子(i)的位置由上述公式决定，通过时间t的变化产生连续的变形和运动，形成类似生物脊骨和触须的视觉效果。具体特点包括：

波动身体：通过参数k和q的复杂组合形成主体波动

触须效果：双曲线项0.3/k在k接近零时产生长触须状结构

螺旋运动：X和Y坐标中的cos(d-t)和sin(d-t)项产生旋转效果

时间演化：所有公式都依赖于时间t，产生连续动画。

四其他数字生命的代码实现

篇幅所限，这里只再给出一个水母的效果实现，其他完整代码可访问gitee直接获取（https://gitee.com/flamingskyline/digital-life）。

4.1 水母

% 水母
figure('Position',[300,50,900,900], 'Color','k');
axes(gcf, 'NextPlot','add', 'Position',[0,0,1,1], 'Color','k');
axis([0, 400, 0, 400])
SHdl = scatter([], [], 1, 'filled','o','w', 'MarkerEdgeColor','none', 'MarkerFaceAlpha',.4);
t = 0;
i = 0:1e4;
x = mod(i, 200);
y = i./43;
k = 5.*cos(x./14).*cos(y./30);
e = y./8 - 13;
d = (k.^2 + e.^2)./59 + 4;
a = atan2(k, e);
while true
    t = t + pi/20;
    q = 60 - 3.*sin(a.*e) + k.*(3 + 4./d.*sin(d.^2 - t.*2));
    c = d./2 + e./99 - t./18;
    SHdl.XData = q.*sin(c) + 200; 
    SHdl.YData = (q + d.*9).*cos(c) + 200;
    drawnow; pause(1e-2)
end