1963：Lorenz 混沌模擬 — 蝴蝶煽動了預報的極限#

副標題：LGP-30 上的一次重跑，跑出了不可預報。

案發現場#

1961 年冬天，MIT 的 Lorenz 想偷懶喝杯咖啡。 Royal McBee LGP-30 正在積分十二變量大氣模式。他把中間結果抄下，捨去小數後幾位重跑。回來一看天氣全變：新軌跡先緊貼，後徹底分岔。現場只有一台每秒幾百次運算的弱電腦。嫌疑人只有一個：被截斷的初值差異。常識說小誤差只致小偏差，預報可無限外推。重跑卻說差之毫釐可謬以千里。 Charney 證明預報可行，Lorenz 發現預報有上限。案件升級：數值天氣的極限在哪裡？

偵查過程#

線索一：從十二維到三維#

Lorenz 把對流簡化到只留三個模態。設 \(x\) 為對流強度， \(y\) 為水平溫差， \(z\) 為垂直廓線畸變。加熱自下而來，耗散與非線性平流互相咬合。參數僅三個：普朗特數 \(\sigma\) 、加熱 \(r\) 、幾何 \(b\) 。 1963 年取 \(\sigma = 10\) ， \(b = 8/3\) ， \(r\) 在 28 附近混沌。變量少到可紙上分析，卻複雜到算不準長期行為。

線索二：截斷初值的分岔實驗#

LGP-30 列印三位小數，記憶體內六位精度。輸入列印值等於給初值動微米級手術。起初差異小於 \(10^{-3}\) ，幾天後開始翻倍。數十天後兩軌跡形同陌路，晴與暴風雨各說各話。截斷更小，分岔只遲到不缺席。誤差近似指數成長，倍增時間是系統內稟屬性。 | 模式天數 | 初值差量級 | 軌跡關係 | |----------|------------|----------| | 0 | 約 \(10^{-3}\) | 幾乎重合 | | 10 | 約 \(10^{-2}\) | 勉強可辨 | | 30 | 量級為一 | 完全無關 |

線索三：三方程的蝴蝶#

第一式是對流與黏性競爭，行內記為 \(dx/dt = \sigma(y-x)\) 。其中 \(x\) 為對流強度， \(y\) 為溫差驅動。第二式是平流與弛豫拉扯，含非線性項 \(xz\) 。第三式是廓線畸變演化，含非線性項 \(xy\) 。完整三方程獨立公式如下。

\[ dx/dt = \sigma(y-x) \]

\[ dy/dt = rx - y - xz \]

\[ dz/dt = xy - bz \]

其中 \(r\) 為加熱強度， \(b\) 為幾何因子， \(z\) 為廓線變量。軌跡永不相交永不重複，形如蝴蝶雙翼。困於有限區域，卻有無窮不重複路徑，即奇異吸引子。

結案報告#

1963 年數值模擬發現決定性混沌，炸毀拉普拉斯妖之夢。數值天氣從此有約兩週預報上限，誤差倍增是物理本質。蝴蝶效應一詞 1972 年定型，成為流行文化符號。氣象界轉向集合預報：報一束可能的未來。 Charney 證明能算，Lorenz 證明不能一直算下去。

證據與工具#

三方程已列上節，此為核心證物，不再重寫。 | 符號 | 意義 | 混沌取值 | |------|------|----------| | \(\sigma\) | 普朗特數 | 取 10 | | \(r\) | 相對瑞利數 | 取 28 | | \(b\) | 幾何因子 | 取 \(8/3\) | | \(x\) | 對流模態 | 蝴蝶一翼座標 | 重現：龍格庫塔積分，步長 0.01 ，初值 \((0, 1, 0)\) 。另一起跑加 \(10^{-5}\) 擾動，驗證差值指數成長。繪三維軌跡確認雙翼，掃 \(r\) 觀察定常分岔混沌轉變。數學後續：李約克定理、費根鮑姆常數、分形維數。因果鏈：Charney 1950 年證明能算，本案證明不能一直算。集合預報由此誕生，一次跑五十個初值，看雲圖分岔。資料同化把觀測揉進模式，盡量延長可預報期。本案檔案編號 CHAOS-1963 ，狀態已結案，蝴蝶仍在飛。 1972 年演講定名蝴蝶效應，巴西蝴蝶掀起德州龍捲風。

補充：程式實作#

咖啡還沒涼，兩條軌跡已經分道揚鑣，本探連追都追不上。

對應程式#

1963-lorenz_butterfly.py

理論呼應#

本文核心是三變量對流模型，混沌取參為 \(r = 28\) ，配上普朗特數與幾何因子。

程式以步長 \(0.01\) 的四階龍格庫塔積分，重演 LGP-30 當年的重跑把戲。

兩條軌跡初值只差 \(1e-8\) ，正是截斷小數重跑的數位翻版。

誤差指數成長、數十天後形同陌路，在此被壓縮成二十個時間單位。

執行方式#

python3 _code/1963-lorenz_butterfly.py

實測輸出#

初值差 1e-8, t=20 分離=6.5109 (驗證 >1: PASS), t=50 分離=11.0736
VERIFY: sep20=6.510889, sep_final=11.073604
ALL CHECKS PASS

關鍵數字有二：二十時間單位後分離為 \(6.5109\) ，遠大於驗證門檻 \(1\) 。

五十時間單位後分離達 \(11.0736\) ，兩條軌跡已分居蝴蝶不同翅膀。

毫釐之差兌成兩位數的分離，這就是預報上限的呈堂證供。

讀者實驗#

把擾動從 \(1e-8\) 改為 \(1e-5\) 或 \(1e-12\) ，或把加熱 \(r\) 從 28 改為 20 ，觀察分離時機與雙翼是否消失。

完整程式如下：

# -*- coding: utf-8 -*-
"""1963 Lorenz 蝴蝶效應 (對應 wiki: 計算模擬學 / 混沌與奇異吸引子)
只用 numpy,固定種子。sigma=10,rho=28,beta=8/3,RK4 dt=0.01 跑 50 時間單位。
兩條初值差 1e-8 的軌跡,驗證 20 時間後分離 > 1 (混沌發散)。
"""
import numpy as np

np.random.seed(0)
sigma, rho, beta = 10.0, 28.0, 8.0 / 3.0
dt = 0.01
T = 50.0
steps = int(T / dt)
split_step = int(20.0 / dt)

def f(s):
    x, y, z = s
    return np.array([sigma * (y - x), x * (rho - z) - y, x * y - beta * z])

def rk4(s):
    k1 = f(s); k2 = f(s + 0.5 * dt * k1)
    k3 = f(s + 0.5 * dt * k2); k4 = f(s + dt * k3)
    return s + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0

s1 = np.array([-8.0, 8.0, 27.0])  # 已在奇異吸引子上,20 時間內即可見指數發散
s2 = s1 + np.array([1e-8, 0.0, 0.0])
sep20 = None
for i in range(steps):
    s1 = rk4(s1); s2 = rk4(s2)
    if i + 1 == split_step:
        sep20 = float(np.linalg.norm(s2 - s1))
sep_final = float(np.linalg.norm(s2 - s1))
ok = sep20 > 1.0
print(f"初值差 1e-8, t=20 分離={sep20:.4f} (驗證 >1: {'PASS' if ok else 'FAIL'}), t=50 分離={sep_final:.4f}")
print(f"VERIFY: sep20={sep20:.6f}, sep_final={sep_final:.6f}")
assert ok, "未觀察到混沌發散"
print("ALL CHECKS PASS")