ローレンツアトラクタで「バタフライ効果」を測ったら、精度を10倍にしても予測できる時間は2.7しか伸びなかった

「蝶の羽ばたきが竜巻を起こす」というたとえで知られるバタフライ効果。この記事では、その核心にある「わずかな初期値の差が、どれくらいの速さで大きくなるか」を数値で確かめる。

この記事では、ローレンツの1963年の論文に登場する、対流を単純化した3変数の方程式を扱う。丸め誤差を発見した当時の気象モデルと、この3変数モデルは区別しておく。

$\frac{dx}{dt} = \sigma(y-x), \quad \frac{dy}{dt} = x(\rho - z) - y, \quad \frac{dz}{dt} = xy - \beta z$

$\sigma=10$、$\rho=28$、$\beta=8/3$。3行しかない。ランダムな要素はどこにも入っていない。

「わずかな差が増幅される」ことは知っている。知りたいのは速さのほうだ。10億分の1のズレは、モデル内の時間でどれくらい経つと、目に見える差になるのか。そして、精度を上げれば予測できる時間はどれだけ伸びるのか。

実験に使ったコードの全文はGitHubに置いています。

まず、アトラクタを見る

4次のルンゲ・クッタ法で解いて、軌道を描く。

2つの渦を行き来しながら、同じ場所を二度は通らない。閉じた輪にならず、かといって発散もせず、有限の領域の中を永遠にさまよう。この形は蝶の羽にも見えるが、ここで測るのは図形の見た目ではなく、初期値の差が増幅される速さだ。

10億分の1のズレが見えるまで

初期値の $x$ を $10^{-9}$ だけずらした2本目を用意して、並走させた。

時間	2本の距離
0	0.000000001
5	0.00000026
10	0.000035
15	0.0025
20	0.027
23.8	1.0

対数目盛でほぼ直線ということは、距離が指数関数的に増えているということだ。$t=23.8$ で距離1.0——アトラクタ全体の大きさが40程度なので、この時点で「同じあたりを走っている2本」ではあるが、もう別の点だ。

動かして見るとこうなる。

前半は完全に重なっていて1本の線にしか見えない。後半で急に分かれる。「じわじわ離れる」ではなく「しばらく何も起きず、あるとき急に」という見え方になるのは、指数関数がそう見えるからだ。

リアプノフ指数を、1回で測ってはいけなかった

この指数的な増大の速さが最大リアプノフ指数 $\lambda$ で、$\rho=28$ のローレンツ系では文献値0.906が知られている。上のグラフの傾きから求めると——0.760。合わない。

原因は測り方だった。1本の摂動ベクトルの伸び方は、アトラクタ上のどこにいるかで大きく変動する。1回の走行だけで傾きを取ると、その区間のたまたまの伸び方を拾ってしまう。

そこで、一定時間ごとに摂動の大きさを元に戻しながら伸び率を積み上げる方法(Benettin法)に変えた。

0.9057。 文献値0.906とほぼ一致した。右のグラフのとおり、最初の100区間では0.889、500区間で0.917、2000区間で0.909と、行ったり来たりしながら収束する。

同じデータ、同じ現象でも、1回の観測から傾きを読むか、何度も測り直して平均するかで0.15違った。前者は0.906に対して16%低い。

精度を上げても、予測できる時間はほとんど伸びない

$\lambda$ が分かると、実務的に一番効く量が出せる。初期誤差 $\delta$ が許容誤差 $\Delta$ に育つまでの時間は

$T = \frac{1}{\lambda}\ln\frac{\Delta}{\delta}$

$\delta$ の対数にしか依存しない。つまり、精度を10倍にしても予測時間は $\ln 10/\lambda$ しか伸びない。実測してみた。

初期のズレ	予測がもつ時間
$10^{-3}$	3.6
$10^{-7}$	18.5
$10^{-11}$	28.5
$10^{-15}$	34.0

精度を10倍にするごとに伸びるのは2.71(長時間平均 $\lambda=0.9057$ を使う近似では $\ln10/\lambda \approx 2.54$。短い区間の傾き0.760を使うと約3.03になる)。$10^{-3}$ から $10^{-15}$ へ、1兆倍精密にして、伸びたのは10倍にも満たない。

ここでの $t=34.0$ は無次元化されたモデル時間で、34時間という意味ではない。また、機械イプシロンは値の尺度によって絶対誤差が変わるため、一律の初期誤差の下限ではない。この実験は初期精度と予測時間の関係を示す例で、実際の天気予報の日数を直接計算したものではない。

アトラクタの次元は、素直に測れた

マンデルブロ集合の境界の次元は理論値2に対して1.09しか出せなかった。ローレンツアトラクタでも測ってみる。今回は箱を数えるのではなく、軌道上の点対のうち距離が $r$ 未満のものの割合 $C(r)$ を使う相関次元という方法にした。

2.050。 文献値の約2.05とほぼ一致した(上の図の左)。両対数プロットが2桁半にわたってきれいな直線になっている。

同じ「フラクタル次元を測る」でも、マンデルブロ境界とは対照的にすんなり出た。アトラクタは軌道が実際に通る場所の集合なので、点を並べればそのまま標本になる。境界のように「1画素より細い枝」を追いかける必要がない。測りやすさは対象の複雑さではなく、標本の取り方で決まっていた。

羽の乗り換えは、一度も予測できなかった

最後に、素朴な疑問を試した。左右どちらの羽を回っているかを0と1で記録すると、記号の列ができる。この列の履歴から「次に乗り換えるか」を当てられるだろうか。

2万点ぶんで乗り換えは3.55%。つまり「乗り換えない」と言い続けるだけで96.45%当たる。

直前 $k$ 個の並びを文脈として、その文脈での乗り換え率が50%を超えたら「乗り換える」と予測する、という単純な予測器を組んだ。$k$ を1から16まで動かした結果は——

どの $k$ でも、乗り換えの再現率は0.0%だった。 50%を超える文脈がひとつも存在しない。予測器は最後まで「乗り換えない」としか言わなかった。

96.45%という正解率だけ見れば優秀に見える。中身は、何も予測していない。

まとめ

• $10^{-9}$ のズレは $t=23.8$ で距離1.0に育つ。増え方は指数的で、前半は完全に重なって見える
• リアプノフ指数は1回の走行から傾きを取ると0.760。何度も測り直して平均すると0.9057(文献値0.906)
• 精度を10倍にしても予測できる時間は2.71しか伸びない。$10^{-15}$ の摂動では、許容距離に達するまでモデル時間で34.0
• 相関次元は2.050(文献値約2.05)。マンデルブロ境界と違い素直に測れた
• 羽の乗り換えは、履歴16個ぶんを使っても再現率0.0%。96.45%の正解率は「乗り換えない」と言い続けた結果

---

効いていたのは、式のどこかに $\ln$ が付いている、それだけだった。

予測できる時間は初期精度の対数に比例して伸びる近似になる。したがって同じ倍率の精度改善なら、同じ時間の延長が目安になる。表の1万倍改善で延長が約15と5.5に分かれたことは、対数則そのものからは説明できない。有限時間の伸び率や丸め誤差などが影響しうるため、この4点だけで原因を確定しない。

仕事で見積もりの精度を上げようとするとき、私は暗黙にこれと逆のことを期待している気がする。調査を2倍丁寧にやれば、見通せる期間も2倍になるだろう、と。実際には、丁寧さは掛け算で効いても、見通せる期間には足し算でしか効かない。

そして羽の乗り換えのほうは、もっと直接的な戒めだった。96.45%当てる予測器を作って、乗り換えは一度も当てていない。知りたかったのは残り3.55%のほうだったのに、正解率という数字はそれを隠す。