七、物理建模
譯者：飛龍
協(xié)議：CC BY-NC-SA 4.0
自豪地采用谷歌翻譯
到目前為止恶导，我們所看到的細胞自動機不是物理模型浆竭；也就是說，他們不打算描述現(xiàn)實世界中的系統(tǒng)惨寿。 但是一些 CA 用作物理模型邦泄。
在本章中，我們考慮一個 CA裂垦，它模擬擴散（散開）并相互反應的化學物質顺囊，這是 Alan Turing 提出的過程，用于解釋一些動物模式如何發(fā)展蕉拢。
我們將試驗一種 CA特碳，它模擬通過多孔材料的滲透液體，例如通過咖啡渣的水晕换。 這個模型是展示相變行為和分形幾何的幾個模型中的第一個午乓，我將解釋這兩者的含義。
本章的代碼位于本書倉庫的chap07.ipynb中届巩。 使用代碼的更多信息硅瞧，請參見第份乒？節(jié)恕汇。
7.1 擴散
1952 年，艾倫圖靈發(fā)表了一篇名為“形態(tài)發(fā)生的化學基礎”的論文或辖，該論文描述了涉及兩種化學物質的系統(tǒng)行為瘾英，它們在空間中擴散并相互反應。 他表明颂暇，這些系統(tǒng)根據(jù)擴散和反應速率產(chǎn)生了廣泛的模式缺谴，并推測像這樣的系統(tǒng)可能是生物生長過程中的重要機制，特別是動物著色模式的發(fā)展耳鸯。
圖靈模型基于微分方程奥洼，但也可以使用細胞自動機來實現(xiàn)鲁沥。
但在我們開始使用圖靈模型之前叁熔，我們先從簡單的事情開始：只有一種化學物質的擴散系統(tǒng)。 我們將使用 2-D CA添谊，其中每個細胞的狀態(tài)是連續(xù)的數(shù)量（通常在 0 和 1 之間），表示化學物質的濃度察迟。
我們將通過比較每個細胞與其鄰居的均值斩狱，來建模擴散過程。 如果中心細胞的濃度超過領域均值扎瓶，則化學物質從中心流向鄰居所踊。 如果中心細胞的濃度較低，則化學物質以另一種方式流動概荷。
以下核計算每個細胞與其鄰居均值之間的差異：
kernel = np.array([[0, 1, 0],
[1,-4, 1],
[0, 1, 0]])
使用np.correlate2d秕岛，我們可以將這個核應用于數(shù)組中的每個細胞：
c = correlate2d(array, kernel, mode='same')
我們將使用一個擴散常數(shù)r，它關聯(lián)了濃度差與流速：
array += r * c
圖 7.1：0乍赫，5 和 10 步后的簡單擴散模型
圖瓣蛀？顯示 CA 的結果，其中n=9, r=0.1雷厂，除了中間的“島”外惋增，初始濃度為 0。 該圖顯示了 CA 的啟動狀態(tài)改鲫，以及 5 步和 10 步之后的狀態(tài)诈皿。 化學物質從中心向外擴散，直到各處濃度相同像棘。
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ù)組代表第一種化學物質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的結果攻礼。乘以ra和rb得出進入或離開每個細胞的擴散速率。
表達式A * B ** 2表示A和B相互反應的比率栗柒。 假設反應消耗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的結果劳闹，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的結果，在A的背景上產(chǎn)生了珊瑚樣的B晦闰。
圖 7.4：1000放祟，2000 和 4000 步之后的反應擴散模型，參數(shù)為f=0.039和k=0.065
圖呻右？顯示了f = 0.039和k = 0.065的結果跪妥。 在類似于有絲分裂的過程中，這些參數(shù)產(chǎn)生的B點生長和分裂声滥，最后形成穩(wěn)定的等距點圖案眉撵。
1952 年以來，觀察和實驗為圖靈猜想提供了一些支持落塑。 目前為止纽疟，看起來許多動物圖案實際上由某種反應擴散過程形成，但尚未證實憾赁。
7.3 滲流
滲流是流體流過半多孔材料的過程污朽。 實例包括巖層中的油，紙中的水和微孔中的氫氣龙考。 滲流模型也用于研究不是嚴格滲濾的系統(tǒng)膘壶，包括流行病和電阻網(wǎng)絡。 請見 http://en.wikipedia.org/wiki/Percolation_theory洲愤。
滲流模型常常用隨機圖來表示颓芭，就像我們在第？章中看到的那樣柬赐，但它們也可以用細胞自動機表示亡问。 在接下來的幾節(jié)中，我們將探索模擬滲流的 2-D CA肛宋。
在這個模型中：
- 最初州藕，每個細胞是概率為
p的“多孔”或者“無孔”，并且除了頂部那行是“濕的”之外酝陈，所有單元都是“干的”床玻。 - 在每個時間步驟中，如果多孔細胞至少有一個濕的鄰居沉帮，它會變濕锈死。 非多孔細胞保持干燥贫堰。
- 模擬運行直至達到不再有細胞改變狀態(tài)的“固定點”。
如果存在從頂部到底部的濕細胞路徑待牵，我們說 CA 具有“滲流簇”其屏。
滲流的一個主要問題是，找到滲流簇的概率以及它如何依賴于p缨该。 這個問題可能會讓你想起第偎行？節(jié)，其中我們計算了隨機 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 中的行數(shù)和列數(shù)。 p是細胞為多孔的概率膨更。
CA 的狀態(tài)存儲在數(shù)組中汗盘，該數(shù)組使用np.random.choice初始化，以概率p選擇 1（多孔）询一，以概率1-p選擇 0（無孔）隐孽。 頂部那行的狀態(tài)設置為 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)設置為 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驮审，并調用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值，這個臨界值似乎是一樣的逝慧。
臨界值附近的行為的快速變化稱為相變昔脯，類似于物理系統(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ù)刻度繪制的結果芋绸。
在每幅圖中媒殉，頂部虛線表示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))
結果是元組列表九火，其中每個元組包含size和size ** 2赚窃，用于比較，以及滲流簇中的細胞數(shù)（不包括頂行中的初始濕細胞）岔激。
圖 7.9：滲流簇中的細胞數(shù)量與 CA 大小
圖勒极？展示了 10 到 100 范圍內的結果。點展示了每個滲流簇中的細胞數(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)下榆俺，森林分形的幾何形狀是什么？ 它的分形維度是多少坞淮？
