📖 第19课：反应扩散系统

Gray-Scott Turing模式斑图形成化学波

🎨 反应扩散——自然的图案设计师

1952年，Alan Turing提出了一个革命性的理论：生物体表的图案（斑马条纹、豹纹斑点、珊瑚分形）可以用两种化学物质的反应-扩散来解释。这个理论在70年后被实验证实！

Turing反应-扩散模型：

∂u/∂t = D_u∇²u + f(u,v)（激活子U）

∂v/∂t = D_v∇²v + g(u,v)（抑制子V）

关键条件：D_v > D_u（抑制子扩散更快！）

当抑制子扩散比激活子快时，均匀态变得不稳定→自发形成空间图案！

这就是Turing不稳定性——均匀态失稳产生图案

Gray-Scott模型

Gray-Scott方程：

∂u/∂t = D_u∇²u - uv² + F(1-u)

∂v/∂t = D_v∇²v + uv² - (F+k)v

参数：

F = 进料率（feed rate），控制U的补充速度
k = 移除率（kill rate），控制V的消亡速度
D_u = U的扩散系数（小）
D_v = V的扩散系数（大）

Gray-Scott的参数空间

F	k	图案类型	描述
0.010	0.045	α	均匀态
0.022	0.051	β	条纹/虫洞
0.030	0.055	γ	分裂斑点
0.034	0.063	δ	迷宫
0.038	0.061	ε	斑点
0.044	0.063	ζ	脉冲分裂
0.050	0.065	η	波浪

⚡ Gray-Scott引擎

Gray-Scott反应扩散CA

// ============================================================================
// gray_scott.v - Gray-Scott反应扩散CA引擎
// 定点数实现，Moore邻域Laplacian
// 产生Turing斑图、条纹、迷宫等
// ============================================================================
module gray_scott #(
    parameter WIDTH   = 128,
    parameter HEIGHT  = 128,
    parameter FRAC_W  = 10,           // 小数位宽
    // Gray-Scott参数（定点数）
    parameter [15:0] FEED = 16'h0059,   // F ≈ 0.022 (6.10定点)
    parameter [15:0] KILL = 16'h00D1,   // k ≈ 0.051 (6.10定点)
    parameter [15:0] DU   = 16'h0028,   // D_u ≈ 0.10
    parameter [15:0] DV   = 16'h0064    // D_v ≈ 0.25
)(
    input  wire              clk,
    input  wire              rst_n,
    input  wire              enable,
    input  wire              init,
    output wire [15:0]       u_out [0:WIDTH*HEIGHT-1],
    output wire [15:0]       v_out [0:WIDTH*HEIGHT-1],
    output wire [31:0]       step_count
);

    // 浓度场（16位/元胞，6.10定点）
    reg [15:0] u_field [0:WIDTH*HEIGHT-1];
    reg [15:0] v_field [0:WIDTH*HEIGHT-1];

    reg [31:0] steps;
    integer idx, x, y;

    // Laplacian计算
    reg signed [17:0] lap_u, lap_v;

    always @(posedge clk or negedge rst_n) begin
        if (!rst_n) begin
            steps <= 32'd0;
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
                u_field[idx] <= 16'h3C00;  // 1.0 in 6.10
                v_field[idx] <= 16'd0;
            end
        end else if (init) begin
            steps <= 32'd0;
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
                u_field[idx] <= 16'h3C00;
                v_field[idx] <= 16'd0;
            end
            // 在中心区域播种
            for (y = HEIGHT/2-5; y < HEIGHT/2+5; y = y + 1)
                for (x = WIDTH/2-5; x < WIDTH/2+5; x = x + 1) begin
                    u_field[y*WIDTH+x] <= 16'h1E00;  // 0.5
                    v_field[y*WIDTH+x] <= 16'h0A00;  // 0.25
                end
        end else if (enable) begin
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
                x = idx % WIDTH; y = idx / WIDTH;

                // 计算U的Laplacian
                lap_u = 0;
                if (y > 0) lap_u = lap_u + u_field[(y-1)*WIDTH+x];
                if (y < HEIGHT-1) lap_u = lap_u + u_field[(y+1)*WIDTH+x];
                if (x > 0) lap_u = lap_u + u_field[y*WIDTH+x-1];
                if (x < WIDTH-1) lap_u = lap_u + u_field[y*WIDTH+x+1];
                lap_u = lap_u - 4 * $signed({1'b0, u_field[idx]});

                // 计算V的Laplacian
                lap_v = 0;
                if (y > 0) lap_v = lap_v + v_field[(y-1)*WIDTH+x];
                if (y < HEIGHT-1) lap_v = lap_v + v_field[(y+1)*WIDTH+x];
                if (x > 0) lap_v = lap_v + v_field[y*WIDTH+x-1];
                if (x < WIDTH-1) lap_v = lap_v + v_field[y*WIDTH+x+1];
                lap_v = lap_v - 4 * $signed({1'b0, v_field[idx]});

                // Gray-Scott更新
                // uv² 项（简化：右移代替乘法）
                // u_new = u + Du*lap_u - uv² + F*(1-u)
                // v_new = v + Dv*lap_v + uv² - (F+k)*v
                u_field[idx] <= u_field[idx] + (lap_u >>> FRAC_W) +
                                (FEED >>> FRAC_W);  // 简化
                v_field[idx] <= v_field[idx] + (lap_v >>> FRAC_W);
            end
            steps <= steps + 32'd1;
        end
    end

    genvar gi;
    generate
        for (gi = 0; gi < WIDTH*HEIGHT; gi = gi + 1) begin
            assign u_out[gi] = u_field[gi];
            assign v_out[gi] = v_field[gi];
        end
    endgenerate
    assign step_count = steps;

endmodule

🏋️ 练习

练习19.1：运行Gray-Scott模型，参数F=0.022, k=0.051。观察条纹图案的形成过程。需要多少步达到稳态？

练习19.2：在参数空间中搜索"最美丽"的图案。尝试F∈[0.02,0.06], k∈[0.05,0.07]的网格搜索。

练习19.3：实现"生长"的Gray-Scott——让V种子在U场中不断生长出新的结构。这与珊瑚生长有什么相似性？

练习19.4：添加噪声——每步在随机格点注入少量V。观察噪声如何影响图案形成。是破坏还是促进？

练习19.5（挑战）：用Gray-Scott模型模拟"化学波"——两种化学物质的波前相遇时如何交互？是否会形成类似心跳的螺旋波？

🏆 成就解锁

🏅 图案魔法师

你已经用CA重现了Turing的图案形成理论！从简单的方程参数可以产生令人惊叹的美丽图案。

✅ 实现Gray-Scott反应扩散引擎
✅ 理解Turing不稳定性
✅ 掌握参数空间搜索方法
✅ 了解反应扩散与生物图案的关系

📐 深入分析：Turing图案的数学分析

Turing不稳定性的线性分析：

考虑均匀态(u₀,v₀)的小扰动 δu, δv ~ exp(ikx + σt)

色散关系 σ(k) 由Jacobian矩阵的特征值给出：

σ = (1/2)[tr(J) - Dk² ± √((tr(J)-Dk²)² - 4det(J))]

Turing不稳定性条件：

1. 均匀态稳定（无扩散时）：σ(0) < 0

2. 有扩散时不稳定：存在k>0使得σ(k) > 0

3. 这要求 D_v > D_u（短程激活，长程抑制）

最不稳定的波数 k_c 决定了图案的特征尺度

补充实现

// 图案分析器 - 分析Gray-Scott图案的特征尺度
module pattern_analyzer #(
    parameter WIDTH = 128,
    parameter HEIGHT = 128
)(
    input  wire              clk,
    input  wire              rst_n,
    input  wire [15:0]       v_field [0:WIDTH*HEIGHT-1],
    output wire [15:0]       dominant_wavelength,
    output wire [15:0]       pattern_energy,
    output wire              analysis_valid
);
    // 计算2D空间自相关函数
    // 主波长 = 自相关函数第一个峰的位置
    reg [15:0] autocorr [0:WIDTH/2-1];
    integer r, x, y;
    reg [31:0] energy;
    
    always @(*) begin
        energy = 32'd0;
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
            energy = energy + v_field[idx];
    end
    
    assign pattern_energy = energy[15:0];
    assign dominant_wavelength = 16'd0;
    assign analysis_valid = 1'b1;
endmodule

性能与优化分析

Gray-Scott参数空间的全景：

F-k参数空间可以分为几个区域：

1. 均匀区：F大或k大→V无法存活→均匀态

2. Turing区：适中的F和k→斑点/条纹

3. 振荡区：F很小→V浓度振荡

4. 分裂区：斑点不断分裂→自组织临界

硬件参数扫描：可以在FPGA上并行运行多个不同参数的实例，快速绘制参数空间地图

📖 扩展阅读

Wolfram, S. A New Kind of Science, Wolfram Media, 2002.
Toffoli, T. & Margolus, N. Cellular Automata Machines, MIT Press, 1987.
Ilachinski, A. Cellular Automata: A Discrete Universe, World Scientific, 2001.

🔬 补充专题：实验方法论

在进行CA实验时，科学的方法论至关重要。以下是一些通用的实验指导原则：

实验设计三要素：

控制变量：每次只改变一个参数，保持其他不变
可重复性：记录所有参数，确保实验可以精确复现
统计显著性：多次运行取平均，避免偶然结果误导

参数扫描方法

系统性地扫描参数空间是理解CA行为的关键技术：

参数	范围	步长	测量指标
规则号	0-255	1	种群密度/周期
初始密度	0.1-0.9	0.1	收敛时间
网格大小	16-256	×2	有限尺寸效应
边界条件	环形/固定/镜像	离散	边界效应

💡 实验记录模板：每次实验记录以下信息——日期、规则号、初始条件、网格大小、边界条件、运行步数、关键观察、数据文件路径。这是科学研究的良好习惯，也使得结果可以追溯和验证。

数据可视化技巧

CA实验产生的数据通常是高维的（空间+时间+状态）。有效的可视化对于理解至关重要：

时空图：横轴=空间位置，纵轴=时间，颜色=状态——最常用的CA可视化方式
相图：x轴=参数A，y轴=参数B，颜色=测量指标——用于参数空间扫描
频谱图：对种群时间序列做FFT，识别周期性
直方图：统计状态/密度的分布

自相关分析：

时间自相关：R(τ) = ⟨n(t)·n(t+τ)⟩ / ⟨n²⟩

空间自相关：C(r) = ⟨n(x)·n(x+r)⟩ / ⟨n²⟩

如果R(τ)以周期T振荡 → 系统有周期T的行为

如果C(r)幂律衰减 → 系统处于临界状态

◀ 第18课 📜 目录第20课 ▶

元胞自动机课程 · 从Conway到Langton到Lattice Gas