📖 第20课：波传播仿真

波动方程弦振动声波干涉衍射

🌊 波——宇宙的通用语言

从水波到声波到光波到引力波，波动现象无处不在。本课用CA模拟波动方程，观察波的传播、反射、干涉和衍射。

波动方程：

∂²u/∂t² = c²∇²u

离散化（中心差分）：

u(t+1) = 2u(t) - u(t-1) + c²Δt²/Δx² × ∇²u(t)

稳定性条件：c·Δt/Δx ≤ 1/√2（2D CFL条件）

⚡ 波动CA引擎

二维波动方程CA

// ============================================================================
// wave_ca.v - 二维波动方程CA引擎
// 双时间层存储，Von Neumann邻域Laplacian
// 支持阻尼、源点和障碍物
// ============================================================================
module wave_ca #(
    parameter WIDTH  = 128,
    parameter HEIGHT = 128,
    parameter FRAC_W = 8,
    parameter [15:0] C_SQ = 16'h0040,    // c² (8.8定点)，约0.25
    parameter [15:0] DAMPING = 16'hFFF8  // 阻尼系数 (8.8定点)，约0.97
)(
    input  wire              clk,
    input  wire              rst_n,
    input  wire              enable,
    input  wire              init,
    // 源点控制
    input  wire [7:0]        src_x,
    input  wire [7:0]        src_y,
    input  wire              src_en,
    input  wire signed [15:0] src_amp,    // 源幅度
    output wire [15:0]       u_curr [0:WIDTH*HEIGHT-1],
    output wire [31:0]       step_count
);

    // 双时间层
    reg signed [15:0] u_prev [0:WIDTH*HEIGHT-1];  // u(t-1)
    reg signed [15:0] u_curr_r [0:WIDTH*HEIGHT-1]; // u(t)
    reg signed [15:0] u_next [0:WIDTH*HEIGHT-1];  // u(t+1)
    reg [31:0] steps;

    integer idx, x, y;
    reg signed [17:0] laplacian;

    always @(*) begin
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
            x = idx % WIDTH; y = idx / WIDTH;

            // Laplacian（Von Neumann邻域，5点模板）
            laplacian = 0;
            if (y > 0)        laplacian = laplacian + u_curr_r[(y-1)*WIDTH+x];
            if (y < HEIGHT-1) laplacian = laplacian + u_curr_r[(y+1)*WIDTH+x];
            if (x > 0)        laplacian = laplacian + u_curr_r[y*WIDTH+x-1];
            if (x < WIDTH-1)  laplacian = laplacian + u_curr_r[y*WIDTH+x+1];
            laplacian = laplacian - 4 * u_curr_r[idx];

            // 波动方程更新
            // u(t+1) = 2*u(t) - u(t-1) + c²*laplacian
            u_next[idx] = 2 * u_curr_r[idx] - u_prev[idx] +
                          (laplacian * $signed(C_SQ) >>> FRAC_W);

            // 阻尼
            u_next[idx] = (u_next[idx] * $signed(DAMPING)) >>> FRAC_W;

            // 源点注入
            if (src_en && x == src_x && y == src_y)
                u_next[idx] = u_next[idx] + src_amp;
        end
    end

    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_prev[idx] <= 16'd0;
                u_curr_r[idx] <= 16'd0;
            end
        end else if (init) begin
            steps <= 32'd0;
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
                u_prev[idx] <= 16'd0;
                u_curr_r[idx] <= 16'd0;
            end
        end else if (enable) begin
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
                u_prev[idx] <= u_curr_r[idx];
                u_curr_r[idx] <= u_next[idx];
            end
            steps <= steps + 32'd1;
        end
    end

    genvar gi;
    generate
        for (gi = 0; gi < WIDTH*HEIGHT; gi = gi + 1)
            assign u_curr[gi] = u_curr_r[gi];
    endgenerate
    assign step_count = steps;

endmodule

🌊 波动现象

现象	设置	观察
圆形波	单点源	同心圆扩散
干涉	两个点源	明暗条纹（建设性/破坏性干涉）
反射	波+固定边界	入射波和反射波叠加
衍射	波通过窄缝	波前弯曲
驻波	两对向传播的波	节点和腹点
多普勒效应	移动源	前方频率高，后方低

🏋️ 练习

练习20.1：模拟单点源产生的圆形波。验证波速c与参数C_SQ的关系。

练习20.2：模拟双缝干涉——在墙上开两个窄缝，让波通过，观察干涉图案。与光学双缝实验对比。

练习20.3：模拟声波在"房间"中的传播——设置四面反射壁和一个声源。观察驻波模式。

练习20.4：实现"透镜"——在某区域增加波速（等效于降低折射率）。观察波前的弯曲。

练习20.5（挑战）：模拟多普勒效应——让声源以一定速度移动，测量不同方向上的表观频率。验证多普勒公式。

🏆 成就解锁

🏅 波动模拟师

物理仿真阶段完成！你已经用CA模拟了流体、扩散、反应和波动——四种最基本的物理过程。

✅ 实现二维波动方程CA
✅ 理解CFL稳定性条件
✅ 观察干涉、衍射、反射
✅ 掌握双时间层更新方法

📐 深入分析：波动方程的能量守恒

波动方程的能量：

总能量 E = (1/2)∫[(∂u/∂t)² + c²|∇u|²]dA

动能：E_k = (1/2)∫(∂u/∂t)²dA

势能：E_p = (1/2)c²∫|∇u|²dA

无阻尼时：E = const（能量守恒）

有阻尼时：dE/dt = -γ∫(∂u/∂t)²dA ≤ 0（能量耗散）

硬件验证：

每步计算E_k和E_p，验证E_k+E_p的守恒性

这是最有效的数值正确性检验——如果能量不守恒，必有实现错误

补充实现

// 能量监测器 - 实时计算波动CA的动能和势能
module energy_monitor #(
    parameter WIDTH = 128,
    parameter HEIGHT = 128,
    parameter FRAC_W = 8
)(
    input  wire              clk,
    input  wire              rst_n,
    input  wire signed [15:0] u_prev [0:WIDTH*HEIGHT-1],
    input  wire signed [15:0] u_curr [0:WIDTH*HEIGHT-1],
    output wire [31:0]       kinetic_energy,
    output wire [31:0]       potential_energy,
    output wire [31:0]       total_energy,
    output wire              energy_conserved
);
    reg [31:0] ek, ep;
    integer idx, x, y;
    reg signed [31:0] du_dt, grad_sq;
    
    always @(*) begin
        ek = 32'd0; ep = 32'd0;
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
            // 动能: (1/2)(u_curr - u_prev)^2
            du_dt = u_curr[idx] - u_prev[idx];
            ek = ek + (du_dt * du_dt) >>> (FRAC_W * 2);
            // 势能: (1/2)c^2|grad u|^2 (简化)
            ep = ep + (u_curr[idx] * u_curr[idx]) >>> (FRAC_W * 2);
        end
    end
    
    assign kinetic_energy = ek;
    assign potential_energy = ep;
    assign total_energy = ek + ep;
    assign energy_conserved = 1'b1;  // 需要与上一步比较
endmodule

性能与优化分析

波动CA的应用：

1. 室内声学：模拟声音在房间中的反射和混响

2. 地震波传播：模拟P波和S波在地球内部的传播

3. 电磁波仿真：FDTD方法的本质就是波动方程的CA化

4. 量子力学：薛定谔方程在形式上类似于扩散+波动方程

📖 扩展阅读

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)幂律衰减 → 系统处于临界状态

🔧 工程实践指南

将CA从理论变为可运行的硬件系统需要解决许多工程细节。以下是基于实践经验的详细指南：

时钟域设计

CA系统通常需要多个时钟域：

clk_ca：CA计算时钟，100-200MHz
clk_vga：VGA像素时钟，25-165MHz
clk_uart：UART波特率时钟，115200Hz
clk_ddr：DDR内存接口时钟，200-400MHz

跨时钟域同步使用双触发器或异步FIFO

资源优化清单

优化项	方法	节省量
邻居计数器	增量更新替代全加法	~40% LUT
状态存储	BRAM替代分布式RAM	~60% FF
规则查找	硬编码XOR替代MUX	~75% LUT（XOR规则）
显示输出	行缓冲替代全帧缓冲	~50% BRAM
边界处理	环形替代固定（零开销）	0

调试技巧

CA系统的调试有其特殊性：

守恒律验证：粒子数守恒、动量守恒、能量守恒——每步检查
已知模式测试：用已知的经典图案（滑翔机、Sierpinski等）验证规则正确性
对偶验证：软件模拟和硬件实现对比，用相同输入产生相同输出
边界测试：特别注意边界条件——环形边界最容易出错
初始条件敏感：不同的初始种子可能暴露不同的bug

⚠️ 常见陷阱：

忘记双缓冲——读写同一状态导致部分元胞使用t+1值
环形边界索引错误——首尾元胞的邻居指向错误位置
定点数溢出——大网格的种群计数器可能溢出
跨时钟域未同步——导致亚稳态和随机错误
VGA时序违反——消隐期内输出数据导致显示异常

性能基准测试

为CA引擎建立性能基准：

关键性能指标：

SPS（Steps Per Second）：每秒完成的CA步数
PPS（Pixel Updates Per Second）：每秒更新的元胞数 = SPS × W × H
延迟：从输入变化到输出响应的时间
功耗：每百万元胞更新所需的能量

理论峰值PPS = f_clk × W × H（全并行）

实际PPS取决于架构——行缓冲约为理论值的1/H

◀ 第19课 📜 目录第21课 ▶

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