1. 布朗粒子运动模拟：从理论到代码的完整实践

布朗运动，这个听起来有点物理课本味道的词，其实离我们并不遥远。从空气中花粉的随机舞动，到金融市场的价格波动，再到微观世界里蛋白质分子的扩散，其背后的核心数学模型都是相通的。作为一名长期和数据、模型打交道的从业者，我经常需要构建这类随机过程的模拟来验证算法、测试系统稳定性，或者直观地理解某个复杂现象。今天，我就来拆解一下“布朗粒子运动模拟”这个项目，它绝不仅仅是生成几条随机轨迹那么简单。我们将从物理图像出发，深入到数值算法的选择、代码实现的细节，再到如何分析模拟结果，并分享一些在实操中容易踩坑的地方。无论你是物理、化学、生物领域的研究者，还是金融工程、机器学习方向的工程师，只要你的工作涉及随机过程或不确定性建模，这篇内容都能给你提供一套可直接复现的“工具箱”。

2. 项目核心：理解布朗运动的物理与数学模型

模拟的第一步是彻底理解你要模拟的对象。布朗运动在物理学中描述的是悬浮在流体中的微小颗粒，由于受到周围流体分子不平衡的碰撞而表现出的无规则运动。在数学上，它被抽象为一种理想的随机过程——维纳过程。

2.1 核心物理图像与数学抽象

布朗粒子的运动有两个关键特征： 无记忆性（马尔可夫性） 和 位移的方差与时间成正比 。无记忆性意味着粒子下一步怎么走，只取决于它现在的位置，跟它过去的历史无关。这直接引出了我们用随机微分方程来描述它的方式。

最核心的方程是朗之万方程。对于在流体中受到粘滞阻力的一维布朗粒子，其运动方程为： m * dv/dt = -γv + η(t) 其中， m 是粒子质量， v 是速度， γ 是阻尼系数（与流体粘度和粒子尺寸有关）， η(t) 是随机力，代表分子碰撞的涨落。

对于大多数模拟场景，特别是当阻尼很大（即粒子运动很快达到稳态）时，我们通常采用 过阻尼近似 ，忽略惯性项 m*dv/dt 。这时，方程简化为： dx/dt = √(2D) * ξ(t) 这里， x 是位置， D 是扩散系数（ D = k_B T / γ ， k_B 是玻尔兹曼常数， T 是温度）， ξ(t) 是高斯白噪声，满足 <ξ(t)> = 0 和 <ξ(t)ξ(t')> = δ(t-t') 。这个方程就是我们将要数值求解的起点。

注意 ：选择过阻尼模型还是惯性模型，取决于你关心的 时间尺度 。如果你关心粒子在微秒甚至更短时间内的瞬态行为（如光学镊子实验），惯性可能很重要。但如果你关心秒、分钟量级的长期扩散行为，过阻尼近似在绝大多数情况下都是足够精确且计算更简单的选择。

2.2 扩散系数D：连接微观与宏观的桥梁

扩散系数 D 是整个模拟中 最重要的物理参数 。它不是一个可以随意设定的数字，而是由系统的物理性质决定的。对于球形粒子在粘性流体中，可以使用斯托克斯-爱因斯坦关系来估算： D = k_B T / (6πηR) 其中， η 是流体动力粘度， R 是粒子半径。

举个例子，在室温（T=300K）的水（η≈0.001 Pa·s）中，一个半径为1微米（1e-6 m）的聚苯乙烯小球，其扩散系数大约为： D = (1.38e-23 * 300) / (6 * 3.14 * 0.001 * 1e-6) ≈ 2.2e-13 m²/s 这个值非常小，意味着粒子扩散得很慢。理解 D 的量级，对于后续设置合理的模拟时间步长和空间尺度至关重要。如果你胡乱设置 D=1 ，模拟出来的粒子可能会像闪电一样飞出去，完全失去物理意义。

3. 数值算法选型与实现细节

有了数学模型，下一步就是如何用计算机来“解”这个随机微分方程。这里的关键是离散化。

3.1 欧拉-丸山法：最常用的入门算法

对于方程 dx/dt = √(2D) * ξ(t) ，最直接、最常用的离散化方法是欧拉-丸山法。在时间步长 Δt 内，粒子位置的更新公式为： x_{n+1} = x_n + √(2D Δt) * W_n 其中， W_n 是一个服从标准正态分布（均值为0，方差为1）的随机数。 √(2D Δt) 这个因子确保了离散化后粒子位移的方差为 2DΔt ，与连续理论的预期一致。

为什么是 √(2D Δt) ？这需要一点推导。从定义出发，在时间 Δt 内，布朗粒子的位移 Δx 应该满足 <Δx> = 0 ， <Δx²> = 2DΔt 。如果我们简单地用 Δx = sqrt(2D) * sqrt(Δt) * W ，那么 <Δx²> = 2D * Δt * <W²> = 2DΔt ，因为 <W²>=1 。这就保证了离散模型在统计性质上与连续模型匹配。

3.2 时间步长Δt的选择：精度与效率的权衡

选择 Δt 是模拟中的第一个“坑”。 Δt 太大，离散误差大，模拟不准确； Δt 太小，计算耗时巨增。

一个实用的经验法则是： Δt 应远小于系统特征时间的平方 。对于自由扩散，一个相关的特征时间是粒子扩散过一个自身半径 R 所需的时间，约为 τ ~ R²/D 。例如，对于上面 R=1μm, D=2.2e-13 m²/s 的粒子， τ ≈ (1e-6)² / (2.2e-13) ≈ 4.5秒 。那么 Δt 选择在0.01秒到0.1秒量级是合理的。

更保险的做法是进行 收敛性测试 ：用不同的 Δt （例如1ms， 10ms， 100ms）运行同一段模拟时间，比较粒子的均方位移（MSD）曲线。如果MSD曲线不再随 Δt 减小而发生显著变化，就说明 Δt 已经足够小。

实操心得 ：在项目初期，我建议先用一个较大的 Δt （比如 τ/10 ）快速验证代码逻辑和整体行为。待核心流程跑通后，再逐步减小 Δt 进行精确模拟和收敛性分析。这样可以避免一开始就陷入漫长的等待，提高调试效率。

3.3 代码实现框架（Python示例）

下面是一个标准的一维布朗运动模拟的核心代码框架。我强烈建议使用 NumPy 进行向量化操作，这比用循环快几个数量级。

import numpy as np
import matplotlib.pyplot as plt

def simulate_brownian_1d(D, delta_t, total_time, x0=0.0, seed=None):
    """
    模拟一维布朗运动。

    参数：
    D: 扩散系数 (单位： length^2 / time)
    delta_t: 时间步长 (单位： time)
    total_time: 总模拟时间 (单位： time)
    x0: 初始位置，默认为0
    seed: 随机数种子，用于复现结果

    返回：
    t: 时间数组
    x: 位置数组
    """
    if seed is not None:
        np.random.seed(seed)

    num_steps = int(total_time / delta_t) + 1
    t = np.linspace(0, total_time, num_steps)

    # 生成随机步长：核心公式
    # 注意：这里生成的是每一步的位移增量
    random_increments = np.random.randn(num_steps - 1) * np.sqrt(2 * D * delta_t)

    # 计算累积位移路径
    x = np.zeros(num_steps)
    x[0] = x0
    x[1:] = x0 + np.cumsum(random_increments)

    return t, x

# 模拟参数设置（基于之前的例子）
D = 2.2e-13  # m^2/s
delta_t = 0.01  # s
total_time = 100.0  # s
x0 = 0.0

# 运行模拟
t, x = simulate_brownian_1d(D, delta_t, total_time, x0, seed=42)

# 可视化一条轨迹
plt.figure(figsize=(10, 5))
plt.plot(t, x, lw=0.5)
plt.xlabel('Time (s)')
plt.ylabel('Position x (m)')
plt.title('Sample 1D Brownian Trajectory')
plt.grid(True, alpha=0.3)
plt.show()

这段代码清晰地体现了核心公式。 np.random.randn 生成标准正态分布随机数，乘以 np.sqrt(2*D*delta_t) 就得到了每一步的随机位移。 np.cumsum 则高效地计算了位置的累积和，即运动轨迹。

4. 从一维到多维及复杂边界处理

一维模拟是基础，但现实世界是三维的，而且粒子往往不是在无限空间中自由运动。

4.1 二维与三维布朗运动模拟

扩展到多维非常简单，因为布朗运动在各个方向上是 独立 的。这意味着我们可以对每个坐标轴（x, y, z）独立地运行一维模拟。

def simulate_brownian_nd(D, delta_t, total_time, dim=2, r0=None, seed=None):
    """
    模拟n维布朗运动。

    参数：
    dim: 维度，2或3
    r0: 初始位置向量，默认为原点
    """
    if seed is not None:
        np.random.seed(seed)
    if r0 is None:
        r0 = np.zeros(dim)

    num_steps = int(total_time / delta_t) + 1
    t = np.linspace(0, total_time, num_steps)

    # 生成(dim, num_steps-1)的随机数组，每一行是一个维度的独立增量
    random_increments = np.random.randn(dim, num_steps - 1) * np.sqrt(2 * D * delta_t)

    # 计算轨迹
    r = np.zeros((dim, num_steps))
    r[:, 0] = r0
    r[:, 1:] = r0[:, np.newaxis] + np.cumsum(random_increments, axis=1)

    return t, r

# 模拟2D布朗运动
t, r = simulate_brownian_nd(D, delta_t, total_time, dim=2, seed=42)
x, y = r[0, :], r[1, :]

plt.figure(figsize=(8, 8))
plt.plot(x, y, lw=0.5, alpha=0.7)
plt.scatter(x[0], y[0], c='green', s=100, label='Start', zorder=5)
plt.scatter(x[-1], y[-1], c='red', s=100, label='End', zorder=5)
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.title('2D Brownian Motion Trajectory')
plt.axis('equal')  # 保证x和y轴比例相同，轨迹不变形
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

axis=1 参数在 np.cumsum 中至关重要，它确保我们在每个维度的时间方向上进行累积，而不是跨维度累积。

4.2 引入边界条件：反射壁与吸收壁

自由空间是理想情况。更多时候，粒子被限制在微腔、通道或细胞膜附近。这就需要处理边界条件。两种最常见的是 反射边界 和 吸收边界 。

• 反射边界 ：粒子碰到边界后像光一样反射回来。实现方法是在每个时间步更新位置后进行检查和修正。
• 吸收边界 ：粒子碰到边界就被“移除”或“捕获”，模拟结束或记录其首次到达时间。

以下是一个在 [0, L] 区间内一维反射边界的实现示例：

def simulate_brownian_1d_reflective(D, delta_t, total_time, L, x0=0.0, seed=None):
    """带反射边界的一维布朗运动模拟。边界在0和L处。"""
    if seed is not None:
        np.random.seed(seed)

    num_steps = int(total_time / delta_t) + 1
    t = np.linspace(0, total_time, num_steps)
    x = np.zeros(num_steps)
    x[0] = x0

    for i in range(1, num_steps):
        # 1. 提议新位置（自由运动）
        dx = np.random.randn() * np.sqrt(2 * D * delta_t)
        x_proposed = x[i-1] + dx

        # 2. 处理反射边界
        # 如果超出右边界L
        while x_proposed > L:
            x_proposed = 2 * L - x_proposed  # 反射
        # 如果超出左边界0
        while x_proposed < 0:
            x_proposed = -x_proposed  # 反射

        x[i] = x_proposed

    return t, x

# 测试反射边界
L = 1e-6  # 1微米宽的区间
t_ref, x_ref = simulate_brownian_1d_reflective(D, 0.001, 10.0, L, x0=L/2, seed=42)
plt.plot(t_ref, x_ref, lw=0.5)
plt.axhline(y=L, color='r', linestyle='--', alpha=0.5, label='Boundary')
plt.axhline(y=0, color='r', linestyle='--', alpha=0.5)
plt.xlabel('Time (s)')
plt.ylabel('Position x (m)')
plt.title('1D Brownian Motion with Reflective Boundaries')
plt.legend()
plt.show()

重要提示 ：注意代码中使用的是 while 循环而非 if 。这是因为一次反射后，粒子可能仍然在边界之外（特别是当 dx 很大时）。使用 while 循环能确保粒子被完全反射回区间内。这是实现反射边界时一个非常经典的细节，用 if 会导致错误的“穿透”边界现象。

5. 模拟结果的分析与验证

生成轨迹只是第一步，从轨迹中提取有物理意义的信息才是关键。 均方位移（Mean Squared Displacement, MSD） 是最核心的分析工具。

5.1 计算均方位移（MSD）

对于一条长度为N的轨迹 {r_0, r_1, ..., r_{N-1}} ，时间间隔为 Δt ，其MSD定义为： MSD(τ = nΔt) = < |r_{i+n} - r_i|² >_i 其中， <...>_i 表示对所有可能的起始点 i 求平均。MSD揭示了粒子扩散的模式。

对于自由扩散，理论上 MSD(τ) = 2d * D * τ ，其中 d 是维度。这是一个 过原点的直线 ，斜率就是 2d*D 。

def calculate_msd(trajectory, max_lag=None):
    """
    计算给定轨迹的均方位移（MSD）。
    轨迹形状：(dim, num_steps) 或 (num_steps,) 对于一维。
    """
    if trajectory.ndim == 1:
        trajectory = trajectory.reshape(1, -1)  # 转换为二维统一处理
    dim, num_steps = trajectory.shape

    if max_lag is None:
        max_lag = num_steps // 4  # 通常取总步长的1/4，以保证统计可靠性

    msd = np.zeros(max_lag)
    time_lags = np.arange(1, max_lag + 1)

    for lag in range(1, max_lag + 1):
        # 计算所有可能的位移差
        disp = trajectory[:, lag:] - trajectory[:, :-lag]
        # 计算平方位移并求平均
        sq_disp = np.sum(disp**2, axis=0)  # 对维度求和，得到每个时间差的平方距离
        msd[lag-1] = np.mean(sq_disp)

    return time_lags, msd

# 对之前模拟的2D轨迹计算MSD
t_lags, msd_vals = calculate_msd(r, max_lag=500)  # r是之前模拟的2D轨迹
# 将时间滞后转换为实际时间
tau = t_lags * delta_t

# 绘制MSD图
plt.figure(figsize=(8, 6))
plt.loglog(tau, msd_vals, 'o-', label='Simulated MSD', markersize=3)
# 绘制理论线
theory_msd = 4 * D * tau  # 2d*D*tau, 其中d=2
plt.loglog(tau, theory_msd, 'r--', label=f'Theory: 4Dt (D={D:.2e})', linewidth=2)
plt.xlabel('Time Lag τ (s)')
plt.ylabel('MSD (m²)')
plt.title('Mean Squared Displacement Analysis')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.show()

如果模拟正确，计算出的MSD数据点应该与红色的理论直线基本重合（在统计涨落范围内）。这是 验证你模拟代码正确性的黄金标准 。

5.2 从MSD中提取扩散系数D

我们可以通过线性拟合MSD-tau曲线的初始线性部分（通常在双对数坐标下看是否斜率为1）来反推扩散系数 D 。

# 选择线性区域进行拟合（通常取前10%-20%的数据点）
linear_region = int(len(tau) * 0.2)
fit_tau = tau[:linear_region]
fit_msd = msd_vals[:linear_region]

# 进行线性拟合：MSD = A * tau， 其中A = 2d*D
A, _ = np.polyfit(fit_tau, fit_msd, 1)  # 一阶多项式拟合
D_fitted = A / 4  # 因为我们是2D模拟，2d=4

print(f"Input Diffusion Coefficient D: {D:.4e} m²/s")
print(f"Fitted Diffusion Coefficient D_fit: {D_fitted:.4e} m²/s")
print(f"Relative Error: {abs(D_fitted - D)/D*100:.2f}%")

如果拟合出的 D_fit 与你输入的 D 相差很大（比如超过10%），你就需要回头检查：时间步长 Δt 是否太大？模拟总时间是否足够长以获得好的统计？MSD计算代码是否有误？

6. 高级主题与性能优化技巧

当模拟大量粒子或复杂环境时，性能和算法稳定性变得至关重要。

6.1 模拟大量粒子：向量化与并行化

如果你需要模拟 N 个独立粒子的运动，千万不要写 for 循环一个个模拟。利用 NumPy 的广播机制进行向量化计算，速度可以快上百倍。

def simulate_many_particles(D, delta_t, total_time, num_particles=1000, dim=2, seed=None):
    """模拟大量独立布朗粒子的轨迹（仅最终位置或全部轨迹）。"""
    if seed is not None:
        np.random.seed(seed)

    num_steps = int(total_time / delta_t) + 1
    # 一次性为所有粒子、所有时间步、所有维度生成随机数
    # 形状：(num_particles, dim, num_steps-1)
    random_shapes = (num_particles, dim, num_steps - 1)
    random_increments = np.random.randn(*random_shapes) * np.sqrt(2 * D * delta_t)

    # 初始位置，假设都在原点
    initial_positions = np.zeros((num_particles, dim, 1))

    # 计算所有粒子的全部轨迹（内存消耗大）
    # trajectories = np.concatenate([initial_positions, np.cumsum(random_increments, axis=2)], axis=2)

    # 更常用的：只计算最终位置（节省内存）
    final_displacements = np.sum(random_increments, axis=2)  # 形状 (num_particles, dim)
    final_positions = final_displacements  # 因为初始位置是0

    # 分析最终位置的分布
    # 例如，计算所有粒子最终位置距离原点的均方距离，理论上应等于 2*dim*D*total_time
    msd_final = np.mean(np.sum(final_positions**2, axis=1))
    print(f"理论最终MSD (2*dim*D*t): {2*dim*D*total_time:.4e}")
    print(f"模拟最终MSD (平均): {msd_final:.4e}")

    return final_positions

# 模拟1000个粒子
final_pos = simulate_many_particles(D, 0.01, 100.0, num_particles=1000, dim=2, seed=42)

对于 num_particles 上万甚至百万的模拟，并且需要完整轨迹时，内存可能成为瓶颈。这时可以考虑：

1. 分块计算 ：将粒子分成若干批，逐批模拟和保存结果。
2. 使用更高效的数据类型 ：如 np.float32 代替默认的 np.float64 ，如果精度允许的话。
3. 并行计算 ：使用 multiprocessing 或 joblib 库将粒子分配给多个CPU核心。但要注意，生成随机数时，每个进程需要有独立的随机种子，否则会导致相关性。

6.2 随机数生成器的选择与种子管理

随机数的质量直接影响模拟结果的可靠性。 NumPy 默认的 np.random.randn 使用的是伪随机数生成器（PRNG）。对于严肃的科学计算，建议使用更现代、统计性质更好的生成器。

# 推荐：使用NumPy的随机数生成器类，便于管理和设置种子
rng = np.random.default_rng(seed=42)  # 使用PCG64算法，比旧版MT19937更好

# 在模拟函数中使用
random_increments = rng.standard_normal((num_steps-1,)) * np.sqrt(2*D*delta_t)

种子管理的重要性 ：设置随机种子（ seed ）可以确保每次运行程序都得到 完全相同 的随机序列，这对于调试和结果复现至关重要。在发布结果或与他人协作时，务必记录下使用的随机种子。

6.3 包含外势场的模拟：朗之万方程的数值解

当粒子处在外力场中（如重力、光学镊子的谐波势阱），运动方程变为： dx/dt = μ * F(x) + √(2D) * ξ(t) 其中 μ 是迁移率（ μ = D / (k_B T) ）， F(x) 是外力。这时，简单的欧拉-丸山法可能不稳定，特别是当势场变化剧烈时。需要采用更稳健的算法，如 欧拉-丸山法（带漂移项） ，但需注意时间步长要非常小。

对于形式为 F(x) = -∇U(x) 的保守力，更新公式为： x_{n+1} = x_n + μ * F(x_n) * Δt + √(2D Δt) * W_n

def simulate_brownian_in_potential(D, delta_t, total_time, potential_force, x0=0.0, kBT=4.1e-21, seed=None):
    """模拟一维势场中的布朗运动。potential_force是一个函数，返回力F(x)。"""
    if seed is not None:
        np.random.seed(seed)
    mu = D / kBT  # 迁移率
    num_steps = int(total_time / delta_t) + 1
    t = np.linspace(0, total_time, num_steps)
    x = np.zeros(num_steps)
    x[0] = x0

    for i in range(num_steps - 1):
        F = potential_force(x[i])
        deterministic_drift = mu * F * delta_t
        random_kick = np.random.randn() * np.sqrt(2 * D * delta_t)
        x[i+1] = x[i] + deterministic_drift + random_kick

    return t, x

# 示例：在简谐势阱 U(x) = 0.5 * k * x^2 中模拟，力 F(x) = -k*x
k = 1e-6  # 势阱刚度，单位 N/m
def harmonic_force(x):
    return -k * x

t_h, x_h = simulate_brownian_in_potential(D, 1e-5, 1.0, harmonic_force, x0=0.0, seed=42)
plt.plot(t_h, x_h, lw=0.5)
plt.xlabel('Time (s)')
plt.ylabel('Position x (m)')
plt.title('Brownian Motion in a Harmonic Potential')
plt.grid(True, alpha=0.3)
plt.show()

警告 ：在势场模拟中，时间步长 Δt 的选择需要格外小心。一个经验法则是： 确定性漂移项 μ*F(x)*Δt 应该远小于特征势场尺度 。例如在谐波势中，特征尺度是势阱宽度 sqrt(k_B T / k) 。通常需要做收敛性测试，确保 Δt 减小到一定程度后，粒子的稳态分布（如位置的概率密度函数）不再变化。

7. 常见问题、调试技巧与实战心得

即使理论清晰，代码简单，在实际模拟中还是会遇到各种问题。下面是我总结的一些典型问题和解决方法。

7.1 模拟结果与理论不符的排查清单

当你发现MSD曲线不是直线，或者拟合出的D值偏差很大时，可以按照以下清单排查：

问题现象	可能原因	检查与解决方法
MSD曲线在短时间（小τ）下偏离理论	时间步长 Δt 太大	进行收敛性测试，逐步减小 Δt ，直到MSD曲线的初始部分与理论线重合。
MSD曲线在长时间（大τ）下波动大，不光滑	总模拟时间太短，统计样本不足	增加模拟总时间 total_time ，或模拟更多条独立轨迹（ num_trajectories ）然后求平均MSD。
MSD曲线整体斜率与理论不符	扩散系数 D 设置错误，或单位混乱	仔细核对 D 的计算公式和单位。确保模拟中 D 、 Δt 、 total_time 单位一致（如全用国际单位制）。
粒子轨迹出现不自然的“跳跃”	随机数生成有问题，或 √(2DΔt) 计算溢出/下溢	检查随机数生成器（用 rng.standard_normal ）。对于极小的 D 或 Δt ， sqrt(2*D*Δt) 可能下溢为0，需使用 np.sqrt(2.0*D*Δt) 确保是浮点数运算。
带边界模拟时，粒子出现在边界外	边界条件实现逻辑错误（如用 if 而不是 while ）	仔细检查边界处理代码，确保粒子被严格限制在域内。对于反射边界，使用 while 循环。
势场模拟中，粒子能量“爆炸式”增长	时间步长 Δt 太大，导致数值不稳定	大幅减小 Δt 。对于 stiff 方程，可能需要使用隐式积分方法（如欧拉-丸山法的隐式变种），但这会复杂很多。

7.2 性能优化实战技巧

1. 向量化是生命线 ：尽可能使用 NumPy 的数组操作代替Python循环。对于 N 个粒子的模拟，生成一个 (N, steps) 的随机数数组，然后用 np.cumsum(axis=1) 一次性计算所有轨迹，比 for 循环快百倍不止。
2. 需要完整轨迹吗？ 很多时候我们只关心统计量（如MSD、分布函数）。对于自由扩散，粒子最终位置的分布是高斯分布，其方差只取决于总时间。这意味着你可以直接根据这个分布抽样得到最终位置，而无需模拟整条路径，节省大量计算和存储。 final_position = np.random.randn(num_particles, dim) * np.sqrt(2*D*total_time) 。
3. 内存与存储 ：模拟长时间、多粒子的轨迹会生成海量数据。考虑只每隔 N 步保存一次位置（下采样），或者实时计算并丢弃轨迹，只保留需要的统计量。使用 np.savez_compressed 保存压缩的数组。
4. 善用随机数种子 ：在调试阶段，固定种子以保证结果可复现。在最终生产运行时，可以使用系统时间或真正的随机源来初始化种子，以确保模拟的独立性。

7.3 从模拟到洞察：如何设计你的模拟实验

布朗运动模拟本身不是目的，它只是工具。在开始写代码前，想清楚这几个问题：

• 科学问题是什么？ 你是想验证扩散理论？研究受限环境对扩散的影响？还是标定实验仪器的参数？
• 需要测量什么量？ MSD？首次通过时间分布？空间概率密度？粒子间的相遇概率？
• 需要多少统计量？ 一条轨迹的随机涨落很大。任何有意义的结论都需要对大量独立轨迹进行系综平均。根据中心极限定理，统计误差通常以 1/√N 衰减， N 是独立样本数（轨迹数）。要想误差减半，你需要4倍的样本。
• 如何与实验对比？ 实验数据有噪声、有漂移、采样频率有限。在模拟中，你需要加入类似的噪声模型、考虑相机的积分时间效应，才能进行公平的比较。

我个人最深刻的一个体会是： 模拟的保真度（时间步长多小、粒子数多少）永远取决于你要回答的问题 。一个用于定性展示布朗运动随机性的动画，用 Δt=0.1s 模拟1个粒子就够了。但若要精确测量在复杂通道中比自由扩散慢15%的扩散系数，你可能需要 Δt=1ms ，模拟上万个粒子，并运行数小时。在动手前做好权衡，能节省大量不必要的计算时间。最后，可视化不仅仅是生成漂亮的图片，更是理解数据和发现错误的重要手段。多画图，多看轨迹、看分布、看随时间的变化，你的直觉会在这个过程中被训练得越来越准。