七、物理建模

到目前为止，我们所看到的细胞自动机不是物理模型；也就是说，他们不打算描述现实世界中的系统。 但是一些 CA 用作物理模型。

在本章中，我们考虑一个 CA，它模拟扩散（散开）并相互反应的化学物质，这是 Alan Turing 提出的过程，用于解释一些动物模式如何发展。

我们将试验一种 CA，它模拟通过多孔材料的渗透液体，例如通过咖啡渣的水。 这个模型是展示相变行为和分形几何的几个模型中的第一个，我将解释这两者的含义。

本章的代码位于本书仓库的chap07.ipynb中。 使用代码的更多信息，请参见第？节。

7.1 扩散

1952 年，艾伦图灵发表了一篇名为“形态发生的化学基础”的论文，该论文描述了涉及两种化学物质的系统行为，它们在空间中扩散并相互反应。 他表明，这些系统根据扩散和反应速率产生了广泛的模式，并推测像这样的系统可能是生物生长过程中的重要机制，特别是动物着色模式的发展。

图灵模型基于微分方程，但也可以使用细胞自动机来实现。

但在我们开始使用图灵模型之前，我们先从简单的事情开始：只有一种化学物质的扩散系统。 我们将使用 2-D CA，其中每个细胞的状态是连续的数量（通常在 0 和 1 之间），表示化学物质的浓度。

我们将通过比较每个细胞与其邻居的均值，来建模扩散过程。 如果中心细胞的浓度超过领域均值，则化学物质从中心流向邻居。 如果中心细胞的浓度较低，则化学物质以另一种方式流动。

以下核计算每个细胞与其邻居均值之间的差异：

kernel = np.array([[0, 1, 0],
                   [1,-4, 1],
                   [0, 1, 0]])

使用np.correlate2d，我们可以将这个核应用于数组中的每个细胞：

c = correlate2d(array, kernel, mode='same')

我们将使用一个扩散常数r，它关联了浓度差与流速：

array += r * c

图 7.1：0，5 和 10 步后的简单扩散模型

图？显示 CA 的结果，其中n=9, r=0.1，除了中间的“岛”外，初始浓度为 0。 该图显示了 CA 的启动状态，以及 5 步和 10 步之后的状态。 化学物质从中心向外扩散，直到各处浓度相同。

7.2 反应扩散

现在我们添加第二种化学物。 我将定义一个新对象ReactionDiffusion，它包含两个数组，每个化学物对应一个：

class ReactionDiffusion(Cell2D):

   def __init__(self, n, m, params):
        self.params = params
        self.array = np.ones((n, m), dtype=float)
        self.array2 = np.zeros((n, m), dtype=float)
        island(self.array2, val=0.1, noise=0.1)

n和m是数组中的行数和列数。 params是参数元组，下面我会解释它。

数组代表第一种化学物质A的浓度，它最初是无处不在的。

array2表示B的浓度，除了中间的一个岛屿，它初始为零，并且由island初始化：

def island(a, val, noise):
    n, m = a.shape
    r = min(n, m) // 20
    a[n//2-r:n//2+r, m//2-r:m//2+r] = val
    a += noise * np.random.random((n, m))

岛的半径r是n或m的二十分之一，以较小者为准。 岛的高度是val，在这个例子中是0.1。 此外，随机均匀噪声（值为 0 到noise）添加到整个数组。

这里是更新数组的step函数：

def step(self):
    """Executes one time step."""
    A = self.array
    B = self.array2
    ra, rb, f, k = self.params

    cA = correlate2d(A, self.kernel, **self.options)
    cB = correlate2d(B, self.kernel, **self.options)

    reaction = A * B**2
    self.array += ra * cA - reaction + f * (1-A)
    self.array2 += rb * cB + reaction - (f+k) * B

参数是

ra：

A的扩散速率（类似于前一节中的r）。

rb：

B的扩散速率。在该模型的大多数版本中，rb约为ra的一半。

f：

进给速率，控制着A添加到系统的速度。

k：

移除速率，控制B从系统中移除的速度。

现在让我们仔细看看更新语句：

reaction = A * B**2
self.array += ra * cA - reaction + f * (1-A)
self.array2 += rb * cB + reaction - (f+k) * B

数组cA和cB是将扩散核应用于A和B的结果。乘以ra和rb得出进入或离开每个细胞的扩散速率。

表达式A * B ** 2表示A和B相互反应的比率。 假设反应消耗A并产生B，我们在第一个方程中减去这个项并在第二个方程中加上它。

表达式f * (1-A)决定A加入系统的速率。 当A接近 0 时，最大进给速率为f。 当A接近 1 时，进给速率下降到零。

最后，表达式(f+k) * B决定B从系统中移除的速率。 当B接近 0 时，该比率变为零。

只要速率参数不太高，A和B的值通常保持在 0 和 1 之间。

图 7.2：1000，2000 和 4000 步之后的反应扩散模型，参数为f=0.035和k=0.057

使用不同的参数，该模型可以产生类似于各种动物身上的条纹和斑点的图案。 在某些情况下，相似性是惊人的，特别是当进给和移除参数在空间上变化时。

对于本节中的所有模拟，ra = 0.5，rb = 0.25。

图？显示了f=0.035和k=0.057的结果，B的浓度以较暗的颜色显示。 有了这些参数，系统就向稳定状态演化，在B的黑色背景上有A的光点。

图 7.3：1000，2000 和 4000 步之后的反应扩散模型，参数为f=0.055和k=0.062

图？显示了f = 0.055和k = 0.062的结果，在A的背景上产生了珊瑚样的B。

图 7.4：1000，2000 和 4000 步之后的反应扩散模型，参数为f=0.039和k=0.065

图？显示了f = 0.039和k = 0.065的结果。 在类似于有丝分裂的过程中，这些参数产生的B点生长和分裂，最后形成稳定的等距点图案。

1952 年以来，观察和实验为图灵猜想提供了一些支持。 目前为止，看起来许多动物图案实际上由某种反应扩散过程形成，但尚未证实。

7.3 渗流

渗流是流体流过半多孔材料的过程。 实例包括岩层中的油，纸中的水和微孔中的氢气。 渗流模型也用于研究不是严格渗滤的系统，包括流行病和电阻网络。 请见 http://en.wikipedia.org/wiki/Percolation_theory。

渗流模型常常用随机图来表示，就像我们在第？章中看到的那样，但它们也可以用细胞自动机表示。 在接下来的几节中，我们将探索模拟渗流的 2-D CA。

在这个模型中：

• 最初，每个细胞是概率为p的“多孔”或者“无孔”，并且除了顶部那行是“湿的”之外，所有单元都是“干的”。
• 在每个时间步骤中，如果多孔细胞至少有一个湿的邻居，它会变湿。 非多孔细胞保持干燥。
• 模拟运行直至达到不再有细胞改变状态的“固定点”。

如果存在从顶部到底部的湿细胞路径，我们说 CA 具有“渗流簇”。

渗流的一个主要问题是，找到渗流簇的概率以及它如何依赖于p。 这个问题可能会让你想起第？节，其中我们计算了随机 ER 图连接的概率。 我们会看到这两个模型之间的几个关系。

我定义了一个新类来表示渗流模型：

class Percolation(Cell2D):

    def __init__(self, n, m, p):
        self.p = p
        self.array = np.random.choice([0, 1], (n, m), p=[1-p, p])
        self.array[0] = 5

n和m是 CA 中的行数和列数。 p是细胞为多孔的概率。

CA 的状态存储在数组中，该数组使用np.random.choice初始化，以概率p选择 1（多孔），以概率1-p选择 0（无孔）。 顶部那行的状态设置为 5，表示一个湿细胞。

在每个时间步骤中，我们使用 4 细胞邻域（不包括对角线）来检查任何多孔细胞是否拥有湿的邻居。 这是核：

kernel = np.array([[0, 1, 0],
                   [1, 0, 1],
                   [0, 1, 0]])

这里是step函数：

correlate2d将邻居的状态相加，如果至少有一个邻居是湿的，那么至少大于 5。 最后一行寻找多孔的细胞，a == 1，并且至少有一个湿邻居，c >= 5，并将它们的状态设置为 5，这代表湿的。

图 7.5：渗流模型的前三个步骤，其中n=10和p=0.5

图？显示了n = 10和p = 0.5的渗流模型的前几个步骤。 非多孔细胞为白色，多孔细胞为浅色，湿细胞为深色。

7.4 相变

现在让我们测试 CA 是否包含渗流簇。

def test_perc(perc):
    num_wet = perc.num_wet()

    num_steps = 0
    while True:
        perc.step()
        num_steps += 1

        if perc.bottom_row_wet():
            return True, num_steps

        new_num_wet = perc.num_wet()
        if new_num_wet == num_wet:
            return False, num_steps

        num_wet = new_num_wet

test_perc接受Percolation对象作为参数。 每次循环中，它都会使 CA 前进一个时间步骤。 它检查底部那行，看看有没有湿的细胞；如果有，它返回True，表示存在渗透簇，以及num_steps，它是到达底部所需的时间步数。

在每个时间步骤中，它还计算湿细胞的数量并检查自上一步以来数量是否增加。 如果没有，我们已经到达了固定点，而没有找到一个渗流簇，所以我们返回False。

为了估计渗流簇的概率，我们生成许多随机初始状态并测试它们：

def estimate_prob_percolating(p=0.5, n=100, iters=100):
    count = 0
    for i in range(iters):
        perc = Percolation(n, p=p)
        flag, _ = test_perc(perc)
        if flag:
            count += 1

    return count / iters

estimate_prob_percolating使用给定的p和n值生成 100 个 CA，并调用test_perc来查看其中有多少个具有渗流簇。 返回值是拥有的 CA 的比例。

当p = 0.55时，渗滤簇的概率接近于 0。p = 0.60时，它约为 70%，而在p = 0.65时，它接近于 1。这种快速转变表明p的临界值接近 0.6。

我们可以更精确地使用随机游走来估计临界值。 从p的初始值开始，我们构造一个Percolation对象并检查它是否具有渗透簇。 如果是这样，p可能太高，所以我们减少它。 如果不是，p可能太低，所以我们增加它。

这里是代码：

def find_critical(p=0.6, n=100, iters=100):
    ps = [p]
    for i in range(iters):
        perc = Percolation(n=n, p=p)
        flag, _ = test_perc(perc)
        if flag:
            p -= 0.005
        else:
            p += 0.005
        ps.append(p)
    return ps

find_critical以p的给定值开始并上下调整，返回值的列表。 当n = 100时，ps的平均值约为 0.59，对于从 50 到 400 的n值，这个临界值似乎是一样的。

临界值附近的行为的快速变化称为相变，类似于物理系统中的相变，例如水在冰点处从液体变为固体的方式。

在处于或接近临界点时，各种各样的系统展示了一组共同的行为和特征。这些行为被统称为临界现象。 在下一节中，我们将探究其中的一个：分形几何。

7.5 分形

为了理解分形，我们必须从维度开始。

对于简单的几何对象，维度根据缩放行为而定义。 例如，如果正方形的边长为l，则其面积为l ** 2。 指数 2 表示正方形是二维的。 同样，如果立方体的边长为l，则其体积为l ** 3，这表示立方体是三维的。

更一般来说，我们可以通过测量一个对象的“尺寸”（通过一些定义），将对象的维度估计为线性度量的函数。

例如，我将通过测量一维细胞自动机的面积（“开”细胞的总数），将它的维度估计为行数的函数。

图 7.6：32 个时间步之后，规则为 20，50 和 18 的一维 CA。

图？展示了三个一维 CA，就像我们在第？节中看到的那样。 规则 20（左）产生一组看似线性的细胞，所以我们预计它是一维的。 规则 50（中）产生类似于三角形的东西，所以我们预计它是二维的。 规则 18（右）也产生类似三角形的东西，但密度不均匀，所以其缩放行为并不明显。

我将用以下函数来估计这些 CA 的维度，该函数计算每个时间步之后的细胞数。 它返回一个元组列表，其中每个元组包含i和i ** 2，用于比较，以及细胞总数。

def count_cells(rule, n=500):
    ca = Cell1D(rule, n)
    ca.start_single()

    res = []
    for i in range(1, n):
        cells = np.sum(ca.array)
        res.append((i, i**2, cells))
        ca.step()

    return res

图 7.7：规则 20，50 和 18 的“开”细胞的数量与时间步数。

图？展示以双对数刻度绘制的结果。

在每幅图中，顶部虚线表示y = i ** 2。 两边取对数，我们得到logy = 2logi。 由于该数字在双对数刻度上，因此直线的斜率为2。

同样，底部的虚线表示y = i。 在双对数刻度上，直线的斜率为 1。

规则 20（左）每两个时间步骤产生三个细胞，所以i步后的细胞总数为y = 1.5 i。 两边取对数，我们得到logy = log1.5 + logi，所以在双对数刻度上，我们期待一条斜率为 1 的线。实际上，线的估计的斜率为 1.01。

规则 50（中）在第i个时间步骤中产生i + 1个新细胞，因此i步之后的细胞总数为y = i ** 2 + i。 如果我们忽略第二项并取两边的对数，我们有logy ~ 2 logi，所以当i变大时，我们预计看到一条斜率为 2 的线。事实上，估计的斜率为 1.97。

最后，对于规则 18（右），估计的斜率大约是 1.57，这显然不是 1，2 或任何其他整数。 这表明规则 18 生成的图案具有“分数维度”；也就是说，它是一个分形。

7.6 分形和渗流模型

图 7.8：p=0.6和n=100, 200, 300的渗流模型

现在让我们回到渗透模型。 图？展示了p = 0.6和n = 100, 200, 300的渗流模型中的湿细胞簇。非正式来说，它们类似于在自然界和数学模型中看到的分形模式。

为了估计它们的分形维度，我们可以运行一系列尺寸的 CA，计算每个渗流簇中湿细胞的数量，然后看看随着我们增加 CA 的大小，细胞计数的规模如何增长。

以下循环运行了模拟：

for size in sizes:
    perc = Percolation(size, p=p)
    flag, _ = test_perc(perc)
    if flag:
        num_filled = perc.num_wet() - size
        res.append((size, size**2, num_filled))

结果是元组列表，其中每个元组包含size和size ** 2，用于比较，以及渗流簇中的细胞数（不包括顶行中的初始湿细胞）。

图 7.9：渗流簇中的细胞数量与 CA 大小

图？展示了 10 到 100 范围内的结果。点展示了每个渗流簇中的细胞数。 拟合这些点的线的斜率大约为 1.85，这表明当p接近临界值时，渗滤簇实际上是分形的。

当p大于临界值时，几乎每个多孔细胞都被填充，因此湿单元的数量仅为p * size ** 2，它的维度为 2。

当p远小于临界值时，湿细胞的数量与 CA 的线性大小成比例，因此它的维度为 1。

7.7 练习

练习 1

在第？节中，我们发现 CA 规则 18 产生了一个分形。 你能找到其他产生分形的一维 CA 吗？

注意：Cell1D.py中的Cell1D对象不会从左边绕到右边，对于某些规则它在边界上创建了手工艺品 [?]。你可能想要使用Wrap1D，它是Cell1D的子类。 它也在Cell1D.py中定义。

练习 2

1990 年，Bak，Chen 和 Tang 提出了一种细胞自动机，它是一种森林火灾的抽象模型。 每个细胞处于三种状态之一：空，被树占用或着火。

CA 的规则是：

• 空细胞以概率p被占用。
• 如果任何一个邻居着火，那么带有树的细胞就会燃烧。
• 即使没有邻居着火，带有树的细胞自发燃烧，概率为f。
• 在下一个时间步骤中，着火的细胞变为空细胞。

编写一个实现这个模型的程序。 你可能想要继承Cell2D。 参数的常用值为p = 0.01和f = 0.001，但你可能想要尝试其他值。

从随机初始条件开始，运行 CA 直到它达到稳定状态，树的数量不再持续增加或减少。

在稳定状态下，森林分形的几何形状是什么？ 它的分形维度是多少？