七、物理建模
原文：Chapter 7 Physical modeling 译者：飞龙 协议：CC BY-NC-SA 4.0 自豪地采用谷歌翻译
到目前为止，我们所看到的细胞自动机不是物理模型；也就是说，他们不打算描述现实世界中的系统。 但是一些 CA 用作物理模型。
在本章中，我们考虑一个 CA，它模拟扩散（散开）并相互反应的化学物质，这是 Alan Turing 提出的过程，用于解释一些动物模式如何发展。
我们将试验一种 CA，它模拟通过多孔材料的渗透液体，例如通过咖啡渣的水。 这个模型是展示相变行为和分形几何的几个模型中的第一个，我将解释这两者的含义。
本章的代码位于本书仓库的chap07.ipynb中。 使用代码的更多信息，请参见第？节。
7.1 扩散
1952 年，艾伦图灵发表了一篇名为“形态发生的化学基础”的论文，该论文描述了涉及两种化学物质的系统行为，它们在空间中扩散并相互反应。 他表明，这些系统根据扩散和反应速率产生了广泛的模式，并推测像这样的系统可能是生物生长过程中的重要机制，特别是动物着色模式的发展。
图灵模型基于微分方程，但也可以使用细胞自动机来实现。
但在我们开始使用图灵模型之前，我们先从简单的事情开始：只有一种化学物质的扩散系统。 我们将使用 2-D CA，其中每个细胞的状态是连续的数量（通常在 0 和 1 之间），表示化学物质的浓度。
我们将通过比较每个细胞与其邻居的均值，来建模扩散过程。 如果中心细胞的浓度超过领域均值，则化学物质从中心流向邻居。 如果中心细胞的浓度较低，则化学物质以另一种方式流动。
以下核计算每个细胞与其邻居均值之间的差异：
代码语言：javascript复制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从系统中移除的速度。
现在让我们仔细看看更新语句：
代码语言：javascript复制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 图连接的概率。 我们会看到这两个模型之间的几个关系。
我定义了一个新类来表示渗流模型：
代码语言：javascript复制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] = 5n和m是 CA 中的行数和列数。 p是细胞为多孔的概率。
CA 的状态存储在数组中，该数组使用np.random.choice初始化，以概率p选择 1（多孔），以概率1-p选择 0（无孔）。 顶部那行的状态设置为 5，表示一个湿细胞。
在每个时间步骤中，我们使用 4 细胞邻域（不包括对角线）来检查任何多孔细胞是否拥有湿的邻居。 这是核：
代码语言：javascript复制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 是否包含渗流簇。
代码语言：javascript复制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_wettest_perc接受Percolation对象作为参数。 每次循环中，它都会使 CA 前进一个时间步骤。 它检查底部那行，看看有没有湿的细胞；如果有，它返回True，表示存在渗透簇，以及num_steps，它是到达底部所需的时间步数。
在每个时间步骤中，它还计算湿细胞的数量并检查自上一步以来数量是否增加。 如果没有，我们已经到达了固定点，而没有找到一个渗流簇，所以我们返回False。
为了估计渗流簇的概率，我们生成许多随机初始状态并测试它们：
代码语言：javascript复制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 / itersestimate_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可能太低，所以我们增加它。
这里是代码：
代码语言：javascript复制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 psfind_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 的大小，细胞计数的规模如何增长。
以下循环运行了模拟：
代码语言：javascript复制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 直到它达到稳定状态，树的数量不再持续增加或减少。
在稳定状态下，森林分形的几何形状是什么？ 它的分形维度是多少？
