七、物理建模
譯者：飛龍
協(xié)議：CC BY-NC-SA 4.0
自豪地采用谷歌翻譯
到目前為止，我們所看到的細胞自動機不是物理模型；也就是說，他們不打算描述現(xiàn)實世界中的系統(tǒng)。 但是一些 CA 用作物理模型。
在本章中，我們考慮一個 CA，它模擬擴散（散開）并相互反應的化學物質(zhì)，這是 Alan Turing 提出的過程，用于解釋一些動物模式如何發(fā)展。
我們將試驗一種 CA，它模擬通過多孔材料的滲透液體，例如通過咖啡渣的水。 這個模型是展示相變行為和分形幾何的幾個模型中的第一個，我將解釋這兩者的含義。
本章的代碼位于本書倉庫的chap07.ipynb中。 使用代碼的更多信息，請參見第？節(jié)。
7.1 擴散
1952 年，艾倫圖靈發(fā)表了一篇名為“形態(tài)發(fā)生的化學基礎(chǔ)”的論文，該論文描述了涉及兩種化學物質(zhì)的系統(tǒng)行為，它們在空間中擴散并相互反應。 他表明，這些系統(tǒng)根據(jù)擴散和反應速率產(chǎn)生了廣泛的模式，并推測像這樣的系統(tǒng)可能是生物生長過程中的重要機制，特別是動物著色模式的發(fā)展。
圖靈模型基于微分方程，但也可以使用細胞自動機來實現(xiàn)。
但在我們開始使用圖靈模型之前，我們先從簡單的事情開始：只有一種化學物質(zhì)的擴散系統(tǒng)。 我們將使用 2-D CA，其中每個細胞的狀態(tài)是連續(xù)的數(shù)量（通常在 0 和 1 之間），表示化學物質(zhì)的濃度。
我們將通過比較每個細胞與其鄰居的均值，來建模擴散過程。 如果中心細胞的濃度超過領(lǐng)域均值，則化學物質(zhì)從中心流向鄰居。 如果中心細胞的濃度較低，則化學物質(zhì)以另一種方式流動。
以下核計算每個細胞與其鄰居均值之間的差異：
kernel = np.array([[0, 1, 0],
[1,-4, 1],
[0, 1, 0]])
使用np.correlate2d，我們可以將這個核應用于數(shù)組中的每個細胞：
c = correlate2d(array, kernel, mode='same')
我們將使用一個擴散常數(shù)r，它關(guān)聯(lián)了濃度差與流速：
array += r * c
圖 7.1：0，5 和 10 步后的簡單擴散模型
圖？顯示 CA 的結(jié)果，其中n=9, r=0.1，除了中間的“島”外，初始濃度為 0。 該圖顯示了 CA 的啟動狀態(tài)，以及 5 步和 10 步之后的狀態(tài)。 化學物質(zhì)從中心向外擴散，直到各處濃度相同。
7.2 反應擴散
現(xiàn)在我們添加第二種化學物。 我將定義一個新對象ReactionDiffusion，它包含兩個數(shù)組，每個化學物對應一個：
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是數(shù)組中的行數(shù)和列數(shù)。 params是參數(shù)元組，下面我會解釋它。
數(shù)組代表第一種化學物質(zhì)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）添加到整個數(shù)組。
這里是更新數(shù)組的step函數(shù)：
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
參數(shù)是
ra：
A的擴散速率（類似于前一節(jié)中的r）。
rb：
B的擴散速率。在該模型的大多數(shù)版本中，rb約為ra的一半。
f：
進給速率，控制著A添加到系統(tǒng)的速度。
k：
移除速率，控制B從系統(tǒng)中移除的速度。
現(xiàn)在讓我們仔細看看更新語句：
reaction = A * B**2
self.array += ra * cA - reaction + f * (1-A)
self.array2 += rb * cB + reaction - (f+k) * B
數(shù)組cA和cB是將擴散核應用于A和B的結(jié)果。乘以ra和rb得出進入或離開每個細胞的擴散速率。
表達式A * B ** 2表示A和B相互反應的比率。 假設(shè)反應消耗A并產(chǎn)生B，我們在第一個方程中減去這個項并在第二個方程中加上它。
表達式f * (1-A)決定A加入系統(tǒng)的速率。 當A接近 0 時，最大進給速率為f。 當A接近 1 時，進給速率下降到零。
最后，表達式(f+k) * B決定B從系統(tǒng)中移除的速率。 當B接近 0 時，該比率變?yōu)榱恪?/p>
只要速率參數(shù)不太高，A和B的值通常保持在 0 和 1 之間。
圖 7.2：1000，2000 和 4000 步之后的反應擴散模型，參數(shù)為f=0.035和k=0.057
使用不同的參數(shù)，該模型可以產(chǎn)生類似于各種動物身上的條紋和斑點的圖案。 在某些情況下，相似性是驚人的，特別是當進給和移除參數(shù)在空間上變化時。
對于本節(jié)中的所有模擬，ra = 0.5，rb = 0.25。
圖？顯示了f=0.035和k=0.057的結(jié)果，B的濃度以較暗的顏色顯示。 有了這些參數(shù)，系統(tǒng)就向穩(wěn)定狀態(tài)演化，在B的黑色背景上有A的光點。
圖 7.3：1000，2000 和 4000 步之后的反應擴散模型，參數(shù)為f=0.055和k=0.062
圖？顯示了f = 0.055和k = 0.062的結(jié)果，在A的背景上產(chǎn)生了珊瑚樣的B。
圖 7.4：1000，2000 和 4000 步之后的反應擴散模型，參數(shù)為f=0.039和k=0.065
圖？顯示了f = 0.039和k = 0.065的結(jié)果。 在類似于有絲分裂的過程中，這些參數(shù)產(chǎn)生的B點生長和分裂，最后形成穩(wěn)定的等距點圖案。
1952 年以來，觀察和實驗為圖靈猜想提供了一些支持。 目前為止，看起來許多動物圖案實際上由某種反應擴散過程形成，但尚未證實。
7.3 滲流
滲流是流體流過半多孔材料的過程。 實例包括巖層中的油，紙中的水和微孔中的氫氣。 滲流模型也用于研究不是嚴格滲濾的系統(tǒng)，包括流行病和電阻網(wǎng)絡(luò)。 請見 http://en.wikipedia.org/wiki/Percolation_theory。
滲流模型常常用隨機圖來表示，就像我們在第？章中看到的那樣，但它們也可以用細胞自動機表示。 在接下來的幾節(jié)中，我們將探索模擬滲流的 2-D CA。
在這個模型中：
- 最初，每個細胞是概率為
p的“多孔”或者“無孔”，并且除了頂部那行是“濕的”之外，所有單元都是“干的”。 - 在每個時間步驟中，如果多孔細胞至少有一個濕的鄰居，它會變濕。 非多孔細胞保持干燥。
- 模擬運行直至達到不再有細胞改變狀態(tài)的“固定點”。
如果存在從頂部到底部的濕細胞路徑，我們說 CA 具有“滲流簇”。
滲流的一個主要問題是，找到滲流簇的概率以及它如何依賴于p。 這個問題可能會讓你想起第？節(jié)，其中我們計算了隨機 ER 圖連接的概率。 我們會看到這兩個模型之間的幾個關(guān)系。
我定義了一個新類來表示滲流模型：
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 中的行數(shù)和列數(shù)。 p是細胞為多孔的概率。
CA 的狀態(tài)存儲在數(shù)組中，該數(shù)組使用np.random.choice初始化，以概率p選擇 1（多孔），以概率1-p選擇 0（無孔）。 頂部那行的狀態(tài)設(shè)置為 5，表示一個濕細胞。
在每個時間步驟中，我們使用 4 細胞鄰域（不包括對角線）來檢查任何多孔細胞是否擁有濕的鄰居。 這是核：
kernel = np.array([[0, 1, 0],
[1, 0, 1],
[0, 1, 0]])
這里是step函數(shù)：
correlate2d將鄰居的狀態(tài)相加，如果至少有一個鄰居是濕的，那么至少大于 5。 最后一行尋找多孔的細胞，a == 1，并且至少有一個濕鄰居，c >= 5，并將它們的狀態(tài)設(shè)置為 5，這代表濕的。
圖 7.5：滲流模型的前三個步驟，其中n=10和p=0.5
圖？顯示了n = 10和p = 0.5的滲流模型的前幾個步驟。 非多孔細胞為白色，多孔細胞為淺色，濕細胞為深色。
7.4 相變
現(xiàn)在讓我們測試 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對象作為參數(shù)。 每次循環(huán)中，它都會使 CA 前進一個時間步驟。 它檢查底部那行，看看有沒有濕的細胞；如果有，它返回True，表示存在滲透簇，以及num_steps，它是到達底部所需的時間步數(shù)。
在每個時間步驟中，它還計算濕細胞的數(shù)量并檢查自上一步以來數(shù)量是否增加。 如果沒有，我們已經(jīng)到達了固定點，而沒有找到一個滲流簇，所以我們返回False。
為了估計滲流簇的概率，我們生成許多隨機初始狀態(tài)并測試它們：
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，并調(diào)用test_perc來查看其中有多少個具有滲流簇。 返回值是擁有的 CA 的比例。
當p = 0.55時，滲濾簇的概率接近于 0。p = 0.60時，它約為 70%，而在p = 0.65時，它接近于 1。這種快速轉(zhuǎn)變表明p的臨界值接近 0.6。
我們可以更精確地使用隨機游走來估計臨界值。 從p的初始值開始，我們構(gòu)造一個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的給定值開始并上下調(diào)整，返回值的列表。 當n = 100時，ps的平均值約為 0.59，對于從 50 到 400 的n值，這個臨界值似乎是一樣的。
臨界值附近的行為的快速變化稱為相變，類似于物理系統(tǒng)中的相變，例如水在冰點處從液體變?yōu)楣腆w的方式。
在處于或接近臨界點時，各種各樣的系統(tǒng)展示了一組共同的行為和特征。這些行為被統(tǒng)稱為臨界現(xiàn)象。 在下一節(jié)中，我們將探究其中的一個：分形幾何。
7.5 分形
為了理解分形，我們必須從維度開始。
對于簡單的幾何對象，維度根據(jù)縮放行為而定義。 例如，如果正方形的邊長為l，則其面積為l ** 2。 指數(shù) 2 表示正方形是二維的。 同樣，如果立方體的邊長為l，則其體積為l ** 3，這表示立方體是三維的。
更一般來說，我們可以通過測量一個對象的“尺寸”（通過一些定義），將對象的維度估計為線性度量的函數(shù)。
例如，我將通過測量一維細胞自動機的面積（“開”細胞的總數(shù)），將它的維度估計為行數(shù)的函數(shù)。
圖 7.6：32 個時間步之后，規(guī)則為 20，50 和 18 的一維 CA。
圖？展示了三個一維 CA，就像我們在第？節(jié)中看到的那樣。 規(guī)則 20（左）產(chǎn)生一組看似線性的細胞，所以我們預計它是一維的。 規(guī)則 50（中）產(chǎn)生類似于三角形的東西，所以我們預計它是二維的。 規(guī)則 18（右）也產(chǎn)生類似三角形的東西，但密度不均勻，所以其縮放行為并不明顯。
我將用以下函數(shù)來估計這些 CA 的維度，該函數(shù)計算每個時間步之后的細胞數(shù)。 它返回一個元組列表，其中每個元組包含i和i ** 2，用于比較，以及細胞總數(shù)。
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：規(guī)則 20，50 和 18 的“開”細胞的數(shù)量與時間步數(shù)。
圖？展示以雙對數(shù)刻度繪制的結(jié)果。
在每幅圖中，頂部虛線表示y = i ** 2。 兩邊取對數(shù)，我們得到logy = 2logi。 由于該數(shù)字在雙對數(shù)刻度上，因此直線的斜率為2。
同樣，底部的虛線表示y = i。 在雙對數(shù)刻度上，直線的斜率為 1。
規(guī)則 20（左）每兩個時間步驟產(chǎn)生三個細胞，所以i步后的細胞總數(shù)為y = 1.5 i。 兩邊取對數(shù)，我們得到logy = log1.5 + logi，所以在雙對數(shù)刻度上，我們期待一條斜率為 1 的線。實際上，線的估計的斜率為 1.01。
規(guī)則 50（中）在第i個時間步驟中產(chǎn)生i + 1個新細胞，因此i步之后的細胞總數(shù)為y = i ** 2 + i。 如果我們忽略第二項并取兩邊的對數(shù)，我們有logy ~ 2 logi，所以當i變大時，我們預計看到一條斜率為 2 的線。事實上，估計的斜率為 1.97。
最后，對于規(guī)則 18（右），估計的斜率大約是 1.57，這顯然不是 1，2 或任何其他整數(shù)。 這表明規(guī)則 18 生成的圖案具有“分數(shù)維度”；也就是說，它是一個分形。
7.6 分形和滲流模型
圖 7.8：p=0.6和n=100, 200, 300的滲流模型
現(xiàn)在讓我們回到滲透模型。 圖？展示了p = 0.6和n = 100, 200, 300的滲流模型中的濕細胞簇。非正式來說，它們類似于在自然界和數(shù)學模型中看到的分形模式。
為了估計它們的分形維度，我們可以運行一系列尺寸的 CA，計算每個滲流簇中濕細胞的數(shù)量，然后看看隨著我們增加 CA 的大小，細胞計數(shù)的規(guī)模如何增長。
以下循環(huán)運行了模擬：
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))
結(jié)果是元組列表，其中每個元組包含size和size ** 2，用于比較，以及滲流簇中的細胞數(shù)（不包括頂行中的初始濕細胞）。
圖 7.9：滲流簇中的細胞數(shù)量與 CA 大小
圖？展示了 10 到 100 范圍內(nèi)的結(jié)果。點展示了每個滲流簇中的細胞數(shù)。 擬合這些點的線的斜率大約為 1.85，這表明當p接近臨界值時，滲濾簇實際上是分形的。
當p大于臨界值時，幾乎每個多孔細胞都被填充，因此濕單元的數(shù)量僅為p * size ** 2，它的維度為 2。
當p遠小于臨界值時，濕細胞的數(shù)量與 CA 的線性大小成比例，因此它的維度為 1。
7.7 練習
練習 1
在第？節(jié)中，我們發(fā)現(xiàn) CA 規(guī)則 18 產(chǎn)生了一個分形。 你能找到其他產(chǎn)生分形的一維 CA 嗎？
注意：Cell1D.py中的Cell1D對象不會從左邊繞到右邊，對于某些規(guī)則它在邊界上創(chuàng)建了手工藝品 [?]。你可能想要使用Wrap1D，它是Cell1D的子類。 它也在Cell1D.py中定義。
練習 2
1990 年，Bak，Chen 和 Tang 提出了一種細胞自動機，它是一種森林火災的抽象模型。 每個細胞處于三種狀態(tài)之一：空，被樹占用或著火。
CA 的規(guī)則是：
- 空細胞以概率
p被占用。 - 如果任何一個鄰居著火，那么帶有樹的細胞就會燃燒。
- 即使沒有鄰居著火，帶有樹的細胞自發(fā)燃燒，概率為
f。 - 在下一個時間步驟中，著火的細胞變?yōu)榭占毎?/li>
編寫一個實現(xiàn)這個模型的程序。 你可能想要繼承Cell2D。 參數(shù)的常用值為p = 0.01和f = 0.001，但你可能想要嘗試其他值。
從隨機初始條件開始，運行 CA 直到它達到穩(wěn)定狀態(tài)，樹的數(shù)量不再持續(xù)增加或減少。
在穩(wěn)定狀態(tài)下，森林分形的幾何形狀是什么？ 它的分形維度是多少？
