「来週の土曜は晴れるでしょうか?」——天気予報は私たちの生活に欠かせませんが、なぜ1週間先になると急に当たらなくなるのでしょうか。スーパーコンピュータが地球大気を細かく刻んで計算しているのに、それでも信頼できるのはせいぜい1〜2週間先まで。この限界は計算機の性能不足ではなく、大気を支配する方程式そのものに刻まれた本質的な性質から生まれます。気象学者エドワード・ローレンツが1963年に発見したその性質こそがカオス(chaos) であり、彼が大気対流の単純化モデルから取り出した3本の微分方程式が、本記事の主役ローレンツ方程式(Lorenz equations) です。
ローレンツ方程式は、変数が $x, y, z$ の3つしかなく、見た目はごく単純です。にもかかわらず、初期値をほんのわずか — 例えば小数点以下6桁目だけ — ずらすだけで、その後の軌道は急速に乖離し、長時間後の振る舞いをまったく予測できなくなります。これが有名な バタフライ効果(butterfly effect):「ブラジルの蝶の羽ばたきがテキサスで竜巻を起こすか?」というローレンツ自身の問いかけが指す現象です。そして驚くべきことに、その軌道を3次元空間に描くと、まさに蝶の羽のような美しい形が浮かび上がります — これがローレンツアトラクタ(Lorenz attractor) です。
カオスは「ロマンチックな数学の好奇心」にとどまりません。現代では、
- カオス工学(chaos engineering) — 通信暗号化、ランダム数生成、レーザーの安定化など、予測不能性を逆手に取った応用
- 神経科学 — 脳波(EEG)の非線形動力学解析、心拍変動のカオス指標による疾患診断
- 気候物理・地球科学 — 気候変動シミュレーションにおける長期挙動の理解、海洋熱塩循環の安定性解析
- 金融工学 — 株価ボラティリティの非線形モデル化、市場のレジーム転換の検出
といった広範な分野で、ローレンツ系で確立された解析手法(分岐解析、リアプノフ指数、ストレンジアトラクタの幾何)が日常的に使われています。
本記事の内容
- ローレンツ方程式がRayleigh-Bénard対流(熱対流)からどのように導出されるかの簡略な道筋
- 3つのパラメータ $\sigma, \rho, \beta$ と固定点の線形安定性解析(ヤコビ行列・固有値)
- $\rho$ の増加にともなうピッチフォーク分岐とサブクリティカル・ホップ分岐
- カオス的振る舞いの定量的指標である初期値鋭敏性と最大リアプノフ指数
scipy.integrate.solve_ivpを用いた高精度RK45積分、3D蝶々型アトラクタ可視化、リアプノフ指数の数値計算- 暗号、神経科学、気候モデリングなど現代のカオス研究への応用展望
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ニュートンの運動法則とは?定義から導出まで丁寧に解説
- 連立常微分方程式と力学系の理論
- 常微分方程式の安定性解析 — 固定点と固有値
- 常微分方程式の数値解法 — オイラー法・ルンゲ-クッタ法
- ルンゲ-クッタ法(RK4)の理論と実装
- リアプノフ安定性と安定性関数
直感 — 蝶々の翼が嵐を呼ぶ
数式に入る前に、ローレンツ系のもつ独特の振る舞いを「絵」でつかんでおきましょう。
水を入れた鍋を弱火で温めると、底の水が暖まって上昇し、表面の水が冷えて沈むという対流ロール(convection roll) が現れます。火を弱いままに保てば、対流は左右どちらかに安定して回り続けます。ところが火力を上げていくと、ある時点でロールが反転するようになります。しばらく右回り、突然左回り、また右回り…… この反転のタイミングは決して周期的ではなく、ほんのわずかな揺らぎで次がいつ起きるかが変わります。これがローレンツ系の本質です。
ローレンツ方程式の解 $(x(t), y(t), z(t))$ を3次元空間に描くと、二枚の翼を持つ蝶のような形が現れます。軌道は片方の翼の周りをしばらくぐるぐる回り、突然反対側の翼に飛び移り、また回り、また飛び移る…… を永遠に繰り返します。重要なのは:
- 軌道は決まった領域(アトラクタ)から離れない — 初期値をどう選んでも、長時間経つと必ずこの蝶々形に吸い寄せられる
- しかし軌道は同じ場所を二度通らない — 周期的にならない
- 初期値を10桁目だけ変えても、しばらく経つと別の翼を飛んでいる — 長期予測は本質的に不可能
これが決定論的カオス(deterministic chaos) です。「決定論的」とは、方程式とパラメータと初期値さえ与えれば未来が完全に決まる、という意味。にもかかわらず実用的には予測不能 — この矛盾めいた性質がローレンツ系を有名にしました。
なぜこのような奇妙なことが起こるのか? その秘密は方程式の非線形項にあります。次節で、ローレンツ方程式がどこから来たのかを見ていきましょう。
ローレンツ方程式の導出 — Rayleigh-Bénard対流からの簡略化
熱対流の方程式
ローレンツ方程式は天下りに与えられたものではなく、水平に広がる薄い流体層(上面が冷たく、下面が温かい)を支配するナビエ-ストークス方程式と熱伝導方程式から、徹底的な簡略化によって取り出された系です。この対流問題をRayleigh-Bénard対流(レイリー-ベナール対流) と呼びます。
出発点は次の3本の偏微分方程式系です(2次元 $(\xi, \zeta)$ で考え、流れ関数 $\psi$ と温度偏差 $\theta$ を用います):
$$ \frac{\partial}{\partial t}\nabla^2 \psi + \frac{\partial(\psi, \nabla^2 \psi)}{\partial(\xi, \zeta)} = \nu \nabla^4 \psi + g\alpha \frac{\partial \theta}{\partial \xi} $$
$$ \frac{\partial \theta}{\partial t} + \frac{\partial(\psi, \theta)}{\partial(\xi, \zeta)} = \frac{\Delta T}{H}\frac{\partial \psi}{\partial \xi} + \kappa \nabla^2 \theta $$
ここで $\nu$ は動粘度、$\kappa$ は熱拡散率、$g$ は重力加速度、$\alpha$ は熱膨張率、$\Delta T$ は上下面の温度差、$H$ は流体層の厚さです。記号 $\partial(f,g)/\partial(\xi,\zeta)$ はヤコビアン(2変数関数の非線形相互作用)を表します。
これらは無限次元の非線形偏微分方程式系で、そのままでは解析的に解けません。ローレンツの天才は、ここで「最も粗いGalerkin近似」を採用したことです。
Galerkin近似 — 3つのモードに限定
流れ関数 $\psi$ と温度偏差 $\theta$ を、対流ロールの最も単純な形を表す三角関数の組で展開します:
$$ \psi(\xi, \zeta, t) = \frac{(1+a^2)\kappa\sqrt{2}}{a}\, X(t)\, \sin\!\Big(\frac{\pi a \xi}{H}\Big) \sin\!\Big(\frac{\pi \zeta}{H}\Big) $$
$$ \theta(\xi, \zeta, t) = \frac{\Delta T R_c}{\pi R}\Big[\sqrt{2}\, Y(t)\, \cos\!\Big(\frac{\pi a \xi}{H}\Big)\sin\!\Big(\frac{\pi \zeta}{H}\Big) – Z(t)\, \sin\!\Big(\frac{2\pi \zeta}{H}\Big)\Big] $$
ここで $a$ は対流ロールのアスペクト比、$R$ はレイリー数(Rayleigh number)、$R_c$ はその臨界値です。展開係数 $X, Y, Z$ はそれぞれ:
- $X(t)$: 対流の流れの強さ(対流ロールの渦度の振幅)
- $Y(t)$: 上昇流と下降流の温度差
- $Z(t)$: 鉛直方向の温度プロファイルの線形性からのずれ(歪み)
を表します。これらを偏微分方程式に代入し、各モードに対する射影方程式を立てると(直交関数の積分は他のモードを消す)、無次元時間 $\tau = (\pi^2/H^2)(1+a^2)\kappa t$ への変数変換を経て、次の3本の常微分方程式が得られます:
$$ \begin{aligned} \frac{dX}{d\tau} &= \sigma(Y – X) \\ \frac{dY}{d\tau} &= \rho X – Y – XZ \\ \frac{dZ}{d\tau} &= XY – \beta Z \end{aligned} $$
これがローレンツ方程式です。3つの無次元パラメータは:
- $\sigma = \nu/\kappa$: プラントル数(Prandtl number)(動粘度と熱拡散率の比)
- $\rho = R/R_c$: 相対レイリー数(熱駆動の強さ)
- $\beta = 4/(1+a^2)$: 対流ロールの形状パラメータ
ローレンツが原論文 Deterministic Nonperiodic Flow (1963) で採用した値は $\sigma = 10$, $\beta = 8/3$, $\rho = 28$ で、これは現在も「標準パラメータ」として定着しています。$\beta = 8/3$ はアスペクト比 $a = 1/\sqrt{2}$ から得られ、これは対流が最も不安定になる「最も増えやすい波数」の選択に対応します。
非線形項の正体
ローレンツ方程式の本質は、$-XZ$ と $+XY$ という2つの非線形項にあります。これらは元の対流方程式におけるヤコビアン項(流れと温度場の移流による相互作用)の名残です。線形項だけ(これらの非線形項を消したもの)であれば解は指数関数の組み合わせとなり、カオスは生まれません。「変数が3つ以上ある」「非線形項を含む」「自律系(時間 $t$ を陽に含まない)である」 — この3条件を満たす最小の系の一つが、ローレンツ方程式なのです。
導出のあらすじが見えたところで、次にこの方程式の振る舞いを定量的に分析していきます。出発点は固定点(fixed point) の探索とその安定性解析です。
パラメータと固定点の安定性解析
固定点を求める
$dX/d\tau = dY/d\tau = dZ/d\tau = 0$ となる点を固定点または平衡点と呼びます。これは時間が経っても動かない点で、力学系の「骨格」を成します。ローレンツ方程式に $0$ を代入すると:
$$ \begin{aligned} \sigma(Y – X) &= 0 \\ \rho X – Y – XZ &= 0 \\ XY – \beta Z &= 0 \end{aligned} $$
第1式から $Y = X$。これを第2式に入れると $X(\rho – 1 – Z) = 0$、第3式から $Z = X^2/\beta$。
ケース1: $X = 0$: $Y = 0$, $Z = 0$。すなわち原点 $\mathbf{P}_0 = (0, 0, 0)$ は常に固定点です。これは対流が完全に止まった状態(熱は伝導のみで運ばれる)に対応します。
ケース2: $X \neq 0$: $\rho – 1 – Z = 0$ より $Z = \rho – 1$。これを $Z = X^2/\beta$ に代入して $X^2 = \beta(\rho – 1)$。よって $\rho > 1$ の場合に限り、追加の固定点が2つ現れます:
$$ \mathbf{P}_{\pm} = \big(\pm\sqrt{\beta(\rho-1)},\ \pm\sqrt{\beta(\rho-1)},\ \rho-1\big) $$
これらは定常的な対流ロール(右回り・左回り)に対応します。$\rho > 1$ で熱駆動が伝導の限界を超えると、対流が立ち上がる — というRayleigh-Bénardの古典的な結果が、固定点解析から直接出てきました。
ヤコビ行列と線形安定性
固定点近傍での挙動を調べるには、ベクトル場 $\mathbf{F}(X, Y, Z) = (\sigma(Y-X), \rho X – Y – XZ, XY – \beta Z)$ のヤコビ行列を計算します:
$$ \mathbf{J}(X, Y, Z) = \begin{pmatrix} -\sigma & \sigma & 0 \\ \rho – Z & -1 & -X \\ Y & X & -\beta \end{pmatrix} $$
線形安定性は、各固定点でのヤコビ行列の固有値の実部の符号で決まります。すべての実部が負なら安定、一つでも正なら不安定です。
原点 $\mathbf{P}_0$ の安定性
$\mathbf{P}_0$ では:
$$ \mathbf{J}(\mathbf{0}) = \begin{pmatrix} -\sigma & \sigma & 0 \\ \rho & -1 & 0 \\ 0 & 0 & -\beta \end{pmatrix} $$
3行目が独立しているので固有値の一つは $\lambda_3 = -\beta < 0$。残りの $2 \times 2$ ブロック
$$ \begin{pmatrix} -\sigma & \sigma \\ \rho & -1 \end{pmatrix} $$
の固有値は特性方程式 $\lambda^2 + (\sigma+1)\lambda – \sigma(\rho-1) = 0$ から:
$$ \lambda_{1,2} = \frac{-(\sigma+1) \pm \sqrt{(\sigma+1)^2 + 4\sigma(\rho-1)}}{2} $$
判別式は常に正(実固有値)。積 $\lambda_1\lambda_2 = -\sigma(\rho-1)$ が負になる条件は $\rho > 1$ で、このとき一方の固有値が正、すなわち $\mathbf{P}_0$ はサドル(不安定)となります。逆に $\rho < 1$ なら2固有値とも負で、原点は安定です。
物理的に解釈すると:$\rho < 1$ では熱駆動が弱すぎて対流が起きず、すべての擾乱は減衰して伝導状態に落ち着く。$\rho = 1$ で対流の閾値を超え、原点が不安定化する — これがRayleigh-Bénardの古典結果です。
対流固定点 $\mathbf{P}_{\pm}$ の安定性
$\mathbf{P}_{\pm} = (\pm c, \pm c, \rho-1)$($c = \sqrt{\beta(\rho-1)}$) でのヤコビ行列の特性方程式は、3次方程式
$$ \lambda^3 + (\sigma + \beta + 1)\lambda^2 + (\rho + \sigma)\beta \lambda + 2\sigma\beta(\rho – 1) = 0 $$
になります。ラウス-フルビッツ判別法を適用して、すべての固有値の実部が負である条件(安定の条件)を求めると:
$$ \rho < \rho_H = \sigma\,\frac{\sigma + \beta + 3}{\sigma - \beta - 1}, \quad (\sigma > \beta + 1) $$
標準パラメータ $\sigma = 10$, $\beta = 8/3$ で計算すると:
$$ \rho_H = 10 \cdot \frac{10 + 8/3 + 3}{10 – 8/3 – 1} = 10 \cdot \frac{47/3}{19/3} = \frac{470}{19} \approx 24.74 $$
つまり、$\rho < 24.74$ では対流固定点 $\mathbf{P}_{\pm}$ が安定で、軌道はそのどちらかに螺旋的に収束します。$\rho > 24.74$ でホップ分岐(Hopf bifurcation) が起き、対流固定点も不安定化 — どの固定点も安定でない状態に突入します。ローレンツが選んだ $\rho = 28$ はこの値をわずかに超えており、まさに「全ての固定点が不安定な領域」でカオス的振る舞いが顕在化します。
固定点とその安定性が見えてきたところで、次に $\rho$ を変化させたときに系全体がどう変わるか — 分岐(bifurcation) の言葉でまとめましょう。
ピッチフォーク分岐とホップ分岐
$\rho = 1$ でのピッチフォーク分岐
$\rho$ をパラメータとして $0$ から徐々に増やしていきます。$\rho < 1$ では原点だけが固定点で安定。$\rho = 1$ を超えた瞬間、原点が不安定化し、同時に2つの新しい固定点 $\mathbf{P}_{\pm}$ が原点から枝分かれするように現れます — これがピッチフォーク分岐(pitchfork bifurcation) です。「ピッチフォーク」(熊手)という名は、$\rho$-$X$ 平面に固定点の位置をプロットしたときに分岐図が三股の熊手に見えることに由来します。
ピッチフォーク分岐は対称性の自発的破れを伴います。元のローレンツ方程式は $(X, Y) \to (-X, -Y)$ の変換で不変(対称)ですが、$\rho > 1$ では系はどちらか一方の対流ロール(右回り or 左回り)を選ばざるを得ません。これは強磁性体が冷えると突然 N極が上か下かを「選ぶ」現象と数学的に同じ構造です。
$\rho = \rho_H \approx 24.74$ でのサブクリティカル・ホップ分岐
$\rho$ をさらに増やすと、対流固定点 $\mathbf{P}_{\pm}$ のまわりの固有値の一つが負の実部から正の実部へと遷移します(複素共役対のまま虚軸を横切る)。これがホップ分岐です。通常のホップ分岐では分岐点で安定な周期軌道(リミットサイクル)が生まれますが、ローレンツ系のホップ分岐はサブクリティカル(subcritical) — 分岐点では不安定なリミットサイクルが消滅します。
サブクリティカル・ホップ分岐の結果、$\rho > \rho_H$ では全ての固定点と周期軌道が不安定となり、軌道は「どこにも落ち着けない」状態になります。それでも軌道は無限大に発散せず、ある有界な領域内に閉じ込められたままです。閉じ込められているのに同じ場所には永遠に戻らない — この複雑な軌道がストレンジアトラクタ(strange attractor) です。
$\rho = 13.926$ でのホモクリニック分岐
実は $\rho_H$ の手前、$\rho \approx 13.926$ でホモクリニック分岐という重要な遷移が起きています。原点 $\mathbf{P}_0$ から出た不安定多様体が、戻ってきて $\mathbf{P}_0$ の安定多様体と接続する(ホモクリニック軌道を形成する)現象です。この分岐を境に、軌道の長時間挙動の質が変わります(過渡的カオスの出現)。
分岐の全体像
標準パラメータ($\sigma = 10$, $\beta = 8/3$)でのローレンツ系の分岐シナリオをまとめると:
| $\rho$ の範囲 | 振る舞い |
|---|---|
| $\rho < 1$ | 原点が大域安定。対流なし |
| $1 < \rho < 13.926$ | $\mathbf{P}_{\pm}$ が安定。軌道は対流ロールへ収束 |
| $13.926 < \rho < 24.06$ | $\mathbf{P}_{\pm}$ は局所安定だが、過渡的カオスが共存 |
| $24.06 < \rho < 24.74$ | $\mathbf{P}_{\pm}$ と非自明アトラクタが共存 |
| $\rho > 24.74$ | ストレンジアトラクタが大域的に支配。決定論的カオス |
| $\rho \gg 28$ | 周期窓やT-point などさらに複雑な構造 |
$\rho = 28$ がいかに絶妙な位置にあるかが見えてきます。ローレンツがこの値を選んだのは偶然ではなく、「カオスが最も鮮明に観察できる領域」であり、現在もカオス研究の標準点として用いられています。
ここまでで「いつカオスが現れるか」がわかりました。次の節では、現れたカオスがどのくらい強いかを定量化する道具 — リアプノフ指数を導入します。
初期値鋭敏性とリアプノフ指数
バタフライ効果の定量化
「初期値をわずかに変えると軌道が大きく食い違う」というカオスの性質を、数式で表現することを考えます。時刻 $0$ で初期値 $\mathbf{x}_0$ と $\mathbf{x}_0 + \delta\mathbf{x}_0$ から出発した2つの軌道の差を $\delta\mathbf{x}(t)$ とすると、カオス系では時間とともに指数的に発散することが経験的に知られています:
$$ |\delta\mathbf{x}(t)| \sim |\delta\mathbf{x}_0|\, e^{\lambda t} $$
この発散率の指数 $\lambda$ を最大リアプノフ指数(largest Lyapunov exponent, LLE) と呼びます。$\lambda > 0$ なら任意に小さい摂動でも有限時間で macroscopic に膨らみ、長期予測が原理的に不可能になります。これがカオスの数学的定義の中核です。
厳密な定義
より厳密には、軌道に沿った接ベクトル $\delta\mathbf{x}(t)$ の時間発展は、ヤコビ行列を用いた線形化方程式
$$ \frac{d}{dt}\delta\mathbf{x} = \mathbf{J}(\mathbf{x}(t))\,\delta\mathbf{x} $$
に従います。これを長時間積分したとき、$\delta\mathbf{x}(t)$ の大きさの増加率を平均したものが最大リアプノフ指数です:
$$ \lambda = \lim_{t \to \infty}\lim_{|\delta\mathbf{x}_0| \to 0}\frac{1}{t}\ln\frac{|\delta\mathbf{x}(t)|}{|\delta\mathbf{x}_0|} $$
ローレンツ系には3つのリアプノフ指数 $\lambda_1 \geq \lambda_2 \geq \lambda_3$ があり、標準パラメータでは:
$$ (\lambda_1, \lambda_2, \lambda_3) \approx (0.906,\ 0,\ -14.572) $$
となることが数値的に知られています。$\lambda_1 > 0$ がカオスの存在、$\lambda_2 = 0$ が軌道方向の中立性(力学系の一般性質)、$\lambda_3 < 0$ がアトラクタの引力を意味します。
予測可能ホライズン
リアプノフ指数から、実用的に予測可能な時間スケールを見積もれます。観測精度を $|\delta\mathbf{x}_0| = \epsilon$、許容誤差を $\Delta$ とすると、誤差が $\Delta$ に到達するまでの時間 $T_{\mathrm{pred}}$ は $\epsilon e^{\lambda T_{\mathrm{pred}}} = \Delta$ から:
$$ T_{\mathrm{pred}} = \frac{1}{\lambda}\ln\frac{\Delta}{\epsilon} $$
ローレンツ系では $\lambda_1 \approx 0.906$ なので、観測精度を1000倍(3桁)向上させても予測可能時間は $\ln(1000)/0.906 \approx 7.6$ 時間単位しか伸びません。地球大気でも事情は同じで、これが「2週間の壁」と呼ばれる気象予報の本質的限界の数学的根拠です。
軌道の収縮と位相空間体積
ローレンツ系には全リアプノフ指数の和に関する興味深い性質があります。ベクトル場の発散
$$ \nabla \cdot \mathbf{F} = -\sigma – 1 – \beta $$
が定数で負(標準パラメータで $-13.667$)。リウヴィルの定理により、位相空間中の体積要素は
$$ V(t) = V(0)\, e^{(\nabla \cdot \mathbf{F})\, t} = V(0)\, e^{-13.667\, t} $$
として時間とともに急速に収縮します。指数の和もこの発散に等しく:
$$ \lambda_1 + \lambda_2 + \lambda_3 = -\sigma – 1 – \beta \approx -13.67 $$
これは「軌道は3次元空間に広がっているにもかかわらず、最終的に体積ゼロのフラクタル集合(ストレンジアトラクタ)に吸い寄せられる」 という事実を表します。実際、ローレンツアトラクタのフラクタル次元(カプラン-ヨーク次元)は
$$ D_{KY} = 2 + \frac{\lambda_1 + \lambda_2}{|\lambda_3|} \approx 2 + \frac{0.906}{14.572} \approx 2.062 $$
と、2次元面より少しだけ厚い「2.06次元」になります。これがストレンジアトラクタの「奇妙さ」の正体です。
理論的な地ならしが整いました。ここからはPythonで実際にローレンツ方程式を積分し、ストレンジアトラクタの可視化と最大リアプノフ指数の数値計算を行ってみましょう。
Python実装 — RK45積分・3Dアトラクタ可視化・リアプノフ指数
scipy.integrate.solve_ivp による高精度積分
冒頭で示した元記事のコードは単純な前進オイラー法でしたが、カオス系では数値誤差そのものがリアプノフ指数で増幅されるため、高精度な積分法を使うのが鉄則です。scipy.integrate.solve_ivp の RK45(陽的ルンゲ-クッタ4-5次、ドルマン-プリンス法)は標準的な選択肢です。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# ローレンツ方程式の右辺
def lorenz_rhs(t, state, sigma=10.0, rho=28.0, beta=8.0/3.0):
x, y, z = state
return [
sigma * (y - x),
x * (rho - z) - y,
x * y - beta * z,
]
# 標準パラメータでの積分
t_span = (0.0, 40.0)
t_eval = np.linspace(*t_span, 10000)
initial_state = [1.0, 1.0, 1.0]
sol = solve_ivp(
lorenz_rhs, t_span, initial_state,
t_eval=t_eval, method="RK45",
rtol=1e-9, atol=1e-12, # カオス系では誤差許容を厳しめに
)
x, y, z = sol.y
print(f"積分成功: {sol.success}, ステップ数: {len(sol.t)}")
print(f"x の範囲: [{x.min():.2f}, {x.max():.2f}]")
print(f"z の範囲: [{z.min():.2f}, {z.max():.2f}]")
このコードは標準パラメータでローレンツ方程式を $t = 0$ から $40$ まで積分します。rtol=1e-9 という厳しい相対誤差設定が重要で、これより緩い設定(例えば既定の 1e-3)では数値誤差が早期に蓄積し、長時間の軌道がアトラクタから外れたり、本来とは違う形状になったりします。出力された $x, z$ の範囲(おおむね $x \in [-20, 20]$, $z \in [0, 50]$)は文献値と一致し、軌道が想定どおりアトラクタに乗っていることが確認できます。
3Dアトラクタの可視化
軌道を3次元空間にプロットして、有名な「蝶々」を描きましょう。色を時間に応じて変化させると、軌道がどう動いていったかが視覚的に追えます。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(11, 9))
ax = fig.add_subplot(111, projection="3d")
# 軌道を時間に応じて色分け
n_pts = len(x)
colors = plt.cm.plasma(np.linspace(0, 1, n_pts - 1))
for i in range(n_pts - 1):
ax.plot(x[i:i+2], y[i:i+2], z[i:i+2], color=colors[i], lw=0.5, alpha=0.85)
# 固定点 P_+, P_- も描画
c = np.sqrt((8.0/3.0) * (28.0 - 1.0))
ax.scatter([c, -c], [c, -c], [27, 27], color="red", s=80,
marker="*", label=r"$\mathbf{P}_{\pm}$ (unstable spirals)", zorder=10)
ax.scatter([0], [0], [0], color="black", s=60, marker="o",
label=r"$\mathbf{P}_0$ (saddle)", zorder=10)
ax.set_xlabel("x"); ax.set_ylabel("y"); ax.set_zlabel("z")
ax.set_title(r"Lorenz attractor ($\sigma=10$, $\rho=28$, $\beta=8/3$)")
ax.legend(loc="upper left")
plt.tight_layout()
plt.savefig("lorenz_attractor.png", dpi=150, bbox_inches="tight")
plt.show()
このプロットから、ローレンツアトラクタの典型的な「蝶」形状が描かれます。軌道は赤い星印で示した2つの不安定スパイラル $\mathbf{P}_{\pm}$ の周りを回りながら、不規則なタイミングで一方の翼から他方へジャンプします。原点(黒丸)はサドル固定点で、軌道は決してそこに留まりません。色のグラデーション(紫→黄)を追うと、軌道がアトラクタ全体を「埋め尽くす」ように分布していることが見て取れます — これがエルゴード性の視覚的表現です。
2つの軌道の発散 — バタフライ効果の実演
初期値をわずかに変えた2軌道を同じ図に描き、いつ・どのように発散するかを確認します。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# 2つの初期値: ほぼ同じだが x 成分のみ 1e-9 だけ違う
state_a = [1.0, 1.0, 1.0]
state_b = [1.0 + 1e-9, 1.0, 1.0]
t_span = (0.0, 35.0)
t_eval = np.linspace(*t_span, 7000)
sol_a = solve_ivp(lorenz_rhs, t_span, state_a, t_eval=t_eval,
method="RK45", rtol=1e-12, atol=1e-14)
sol_b = solve_ivp(lorenz_rhs, t_span, state_b, t_eval=t_eval,
method="RK45", rtol=1e-12, atol=1e-14)
# 軌道の差(ユークリッド距離)
diff = np.sqrt(np.sum((sol_a.y - sol_b.y) ** 2, axis=0))
fig, axes = plt.subplots(2, 1, figsize=(11, 8), sharex=True)
axes[0].plot(sol_a.t, sol_a.y[0], label="trajectory A: x0=1.0", lw=0.9)
axes[0].plot(sol_b.t, sol_b.y[0], label="trajectory B: x0=1.0+1e-9", lw=0.9, alpha=0.8)
axes[0].set_ylabel("x(t)")
axes[0].legend(); axes[0].grid(alpha=0.3)
axes[0].set_title("Two trajectories with initial separation $10^{-9}$")
axes[1].semilogy(sol_a.t, diff, color="crimson", lw=1.2,
label=r"$|\delta\mathbf{x}(t)|$")
# 理論的な指数発散ライン (lambda ~ 0.906)
axes[1].semilogy(sol_a.t, 1e-9 * np.exp(0.906 * sol_a.t), "k--",
lw=1.2, label=r"$10^{-9} e^{0.906 t}$ (theory)")
axes[1].set_xlabel("t"); axes[1].set_ylabel("Euclidean distance")
axes[1].legend(); axes[1].grid(alpha=0.3, which="both")
axes[1].set_ylim(1e-12, 1e2)
plt.tight_layout()
plt.savefig("lorenz_divergence.png", dpi=150, bbox_inches="tight")
plt.show()
上段の $x(t)$ の時系列を見ると、$t \lesssim 10$ までは2つの軌道はほぼ完全に重なって見えますが、$t \approx 15$ あたりから明確にずれ始め、それ以降は完全に独立した時系列として振る舞います。下段の対数プロットがその発散をクリアに示しています — 赤い実線(軌道差)は、初期の $10^{-9}$ から 理論的な指数増加曲線 $10^{-9}e^{0.906t}$(黒破線)とほぼ平行に増加 し、約 $t \approx 25$ でアトラクタのサイズ($O(10)$)に飽和します。
この「最初は同じだが、ある時刻から急に違う」という挙動こそがバタフライ効果の数値的実体です。蝶の羽ばたきが大気の状態をわずかに揺らがせた瞬間、その揺らぎは大気のローレンツ的力学によって毎秒 $e^{0.906}$ 倍(約2.5倍)のオーダーで増幅され、数日後には全く別の天気を生み出す可能性がある — そう解釈できます。
最大リアプノフ指数の数値計算 — Bennettin 法
最大リアプノフ指数を数値的に計算する標準的な手法がBennettinアルゴリズム(または「再正規化付き軌道追跡法」)です。基本軌道の傍に微小摂動軌道を置き、(1)両方を短時間積分、(2)摂動の大きさを測る、(3)摂動を初期サイズに再正規化、を繰り返します。
import numpy as np
from scipy.integrate import solve_ivp
def compute_lyapunov(rhs, x0, t_total=200.0, dt_renorm=0.5, eps=1e-8,
transient=50.0, rtol=1e-9, atol=1e-12):
"""最大リアプノフ指数の Bennettin 法による推定。
rhs: ODE 右辺 f(t, x)
x0: 基準軌道の初期値
dt_renorm: 再正規化までの時間間隔
eps: 摂動の固定サイズ
transient: アトラクタに乗るまで捨てる時間
"""
# まず transient まで進めて、アトラクタ上の出発点を取得
sol = solve_ivp(rhs, (0, transient), x0, method="RK45",
rtol=rtol, atol=atol, t_eval=[transient])
x_ref = sol.y[:, -1].copy()
# 基準軌道から eps だけずれた摂動軌道の初期点
x_pert = x_ref + np.array([eps, 0.0, 0.0])
n_steps = int(t_total / dt_renorm)
log_growth_sum = 0.0 # 対数増加率の累計
for k in range(n_steps):
# 両軌道を dt_renorm だけ進める
sol_ref = solve_ivp(rhs, (0, dt_renorm), x_ref, method="RK45",
rtol=rtol, atol=atol, t_eval=[dt_renorm])
sol_pert = solve_ivp(rhs, (0, dt_renorm), x_pert, method="RK45",
rtol=rtol, atol=atol, t_eval=[dt_renorm])
x_ref = sol_ref.y[:, -1]
x_pert = sol_pert.y[:, -1]
# 距離を測り、対数比を累積、それから eps に再正規化
d = np.linalg.norm(x_pert - x_ref)
log_growth_sum += np.log(d / eps)
x_pert = x_ref + (x_pert - x_ref) * (eps / d)
return log_growth_sum / (n_steps * dt_renorm)
lle = compute_lyapunov(lorenz_rhs, [1.0, 1.0, 1.0],
t_total=200.0, dt_renorm=0.5, eps=1e-8)
print(f"推定された最大リアプノフ指数: lambda ~ {lle:.4f}")
print(f"文献値: lambda ~ 0.9056")
実行すると lambda ~ 0.90 付近の値が出力されます。文献値 $0.9056$ と小数2桁まで一致しており、Bennettin法が正しく機能していることがわかります。なぜこの手法が正しいかというと、各ステップで $d/\epsilon = e^{\lambda \Delta t}$ となるはずなので、$\ln(d/\epsilon) = \lambda \Delta t$ を全ステップで平均すれば $\lambda$ が抽出されるからです。再正規化を行わないと摂動がすぐに飽和してしまい、$\ln$ の中身が頭打ちになるため、定期的なリセットが本質的に重要です。
パラメータ $\rho$ を変えた分岐ダイアグラム
最後に、$\rho$ を $0$ から $30$ まで変化させたときに、軌道のピーク値 $z_{\max}$ がどう変わるかを描くポアンカレ分岐図を作ります。これは分岐シナリオを実験的に確認する古典的な手法です。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def collect_peaks(rho, n_peaks=80, t_skip=80.0, t_max=300.0):
"""各 rho について、x(t) の局所最大値を集める"""
sol = solve_ivp(
lambda t, s: lorenz_rhs(t, s, rho=rho),
(0.0, t_max), [1.0, 1.0, 1.0],
method="RK45", rtol=1e-9, atol=1e-12,
dense_output=True,
)
# transient を捨ててから dense で細かくサンプリング
t_dense = np.linspace(t_skip, t_max, 30000)
x_dense = sol.sol(t_dense)[0]
# 局所最大値を検出
peaks = []
for i in range(1, len(x_dense) - 1):
if x_dense[i] > x_dense[i-1] and x_dense[i] > x_dense[i+1]:
peaks.append(x_dense[i])
if len(peaks) >= n_peaks:
break
return np.array(peaks)
rhos = np.linspace(0.5, 30.0, 200)
fig, ax = plt.subplots(figsize=(11, 6))
for rho in rhos:
peaks = collect_peaks(rho)
if len(peaks) > 0:
ax.plot(np.full_like(peaks, rho), peaks, "k.", markersize=0.6, alpha=0.5)
# 重要な分岐点を縦線で
for rho_b, label in [(1.0, r"$\rho=1$ pitchfork"),
(13.926, r"$\rho \approx 13.93$ homoclinic"),
(24.74, r"$\rho \approx 24.74$ Hopf")]:
ax.axvline(rho_b, color="red", ls="--", alpha=0.6, lw=1.0)
ax.text(rho_b, ax.get_ylim()[1] * 0.95, label, rotation=90,
color="red", fontsize=9, va="top", ha="right")
ax.set_xlabel(r"$\rho$ (Rayleigh parameter)")
ax.set_ylabel(r"local maxima of $x(t)$")
ax.set_title("Bifurcation diagram of the Lorenz system")
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("lorenz_bifurcation.png", dpi=150, bbox_inches="tight")
plt.show()
この分岐図から、ローレンツ系の遷移シナリオが一目で読み取れます。$\rho < 1$ では $x$ のピークは $0$(原点に収束)。$\rho \in (1, 13.93)$ では2本の細い線(対流固定点 $\mathbf{P}_{\pm}$ への螺旋的収束のピーク)。$\rho \in (13.93, 24.74)$ では細線とカオス的散布が共存する複雑な領域。$\rho > 24.74$ ではピーク値が広く散布(カオス本格化)し、所々に周期窓(細い帯)が混在します。理論で議論した分岐点が縦の赤破線とぴたりと整合していることが視覚的に確認でき、力学系理論と数値実験の対応関係がはっきり見えます。
異なる $\rho$ でのアトラクタ形状の比較
最後に、$\rho$ の値を変えたときにアトラクタの形がどう変わるかを並べて表示します。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(14, 4.5))
rho_values = [15.0, 24.0, 28.0, 100.0]
titles = [r"$\rho=15$ (固定点)", r"$\rho=24$ (境界)",
r"$\rho=28$ (カオス)", r"$\rho=100$ (周期窓)"]
for idx, (rho, title) in enumerate(zip(rho_values, titles), 1):
sol = solve_ivp(lambda t, s: lorenz_rhs(t, s, rho=rho),
(0, 50), [1.0, 1.0, 1.0], method="RK45",
rtol=1e-9, atol=1e-12, t_eval=np.linspace(20, 50, 5000))
ax = fig.add_subplot(1, 4, idx, projection="3d")
ax.plot(sol.y[0], sol.y[1], sol.y[2], lw=0.4, color="steelblue")
ax.set_title(title); ax.set_xlabel("x"); ax.set_ylabel("y"); ax.set_zlabel("z")
plt.tight_layout()
plt.savefig("lorenz_rho_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
このコードを実行すると、$\rho$ の値ごとに全く異なる軌道形状が並びます。$\rho = 15$ では軌道は片方の不動点に巻きついた小さなスパイラル、$\rho = 24$ では対流固定点を回りつつ稀に乗り換える過渡的カオス、$\rho = 28$ では有名な左右対称の蝶々形ストレンジアトラクタ、$\rho = 100$ では再び規則的な周期軌道(周期窓内)が現れます。「カオスは $\rho$ を上げる一方ではなく、上げたり下げたりすると消えたり現れたりする」という分岐構造の複雑さが、絵として実感できます。
ローレンツ系の応用と現代カオス研究
ローレンツ方程式は1963年の提唱から60年以上経った今もなお、カオス研究の「標準モデル」として使われ続けています。その応用範囲は当初の気象学を遥かに超えて広がっています。
気候・大気科学への応用: ローレンツ系は単純化された大気モデルとして、現在も気候モデルの長期挙動の理解に使われます。アンサンブル予報(複数の初期値からの予報を統計平均する手法)は、初期値鋭敏性に対する実用的な対処法として、ECMWF や気象庁の現業予報で標準化されています。ローレンツ系で確立された「軌道集団の発散度合いから予測可能時間を見積もる」という考え方が、その理論的基盤を成しています。
暗号・乱数生成: $\rho > 24.74$ でのローレンツ系はエルゴード的かつ初期値鋭敏で、軌道はアトラクタを「広く」訪れます。これを利用したカオス暗号は、初期値とパラメータを鍵として擬似乱数列を生成する暗号方式として研究されています。同期化されたローレンツ系を送受信に用いるカオス同期通信も提案されており、傍受者が同じパラメータを知らない限り信号を復号できない、というアイデアです。
神経科学: 脳波(EEG)や心拍変動(HRV)の時系列は非線形動力学的構造を持ち、リアプノフ指数や相関次元といったカオス指標で解析されます。てんかん発作の予兆検出、麻酔深度のモニタリング、統合失調症や鬱病のバイオマーカーとして、これらの指標は実臨床に近づいています。ローレンツ系で開発された Bennettin法は、これらの応用で標準ツールキットの一部です。
化学反応: Belousov-Zhabotinsky反応(BZ反応)などの振動化学反応は、ローレンツに似た低次元カオスを示します。反応器内の温度・濃度の振動パターンから、化学工学的に安定な運転条件を設計する際の指針になります。
レーザー物理: 半導体レーザーや CO₂ レーザーの強度ダイナミクスは、ローレンツ方程式の構造と数学的に等価であることが知られています(Haken-Lorenz モデル)。レーザーが安定発振するか、カオス的に振動するかは、ローレンツ系の分岐シナリオで予測できます。光ファイバ通信のジッタ抑制や 超短パルスレーザの設計で活用されています。
人工知能との融合: 近年は Physics-Informed Neural Networks (PINNs) や Neural ODE を用いてカオス系を学習・予測する研究が活発です。ニューラルネットがリアプノフ指数を保ったまま長時間予測できるかは未解決の挑戦的問題で、気候モデルの代理学習にも応用が広がっています。ローレンツ系はこれらの手法のベンチマーク問題として標準的に用いられます。
機械学習のリザバー計算: ローレンツ系の出力を学習データとして、リザバー(乱雑に結合された再帰ニューラルネット)が将来を予測する Reservoir Computing は、$1/\lambda_1$ 程度の時間スケールまでは驚くほど正確な予測ができることが報告されています。これは「カオス系であってもリアプノフ時間内なら学習可能」という新しい知見をもたらしました。
金融・経済: 株価や為替の時系列に決定論的カオスが含まれるかは古典的な論争ですが、少なくとも非線形性は確実に存在し、ローレンツ系で発展した手法がボラティリティ予測やレジーム検出に応用されています。
これらの応用の根底にあるのは、「自由度がたった3つの単純なモデルで、複雑なシステムの本質的振る舞いを捉えられる」というローレンツのもう一つの教訓です。高次元の現実問題を低次元のローレンツ系に写像し、そこで得た洞察を元の問題に持ち帰る — この「最小モデル」のアプローチは、現代の応用数学に深く根付いています。
まとめ
本記事では、気象モデルから生まれたローレンツ方程式とローレンツアトラクタについて、導出・分岐・初期値鋭敏性・リアプノフ指数までを通して解説しました。
- ローレンツ方程式の起源: Rayleigh-Bénard対流の偏微分方程式を3モードに Galerkin近似することで $\dot X = \sigma(Y-X), \dot Y = \rho X – Y – XZ, \dot Z = XY – \beta Z$ が得られた。
- 3つのパラメータ: $\sigma$(プラントル数), $\rho$(相対レイリー数), $\beta$(対流形状)。標準値は $\sigma=10, \beta=8/3, \rho=28$。
- 固定点と安定性: 原点 $\mathbf{P}_0$ と対流固定点 $\mathbf{P}_{\pm} = (\pm\sqrt{\beta(\rho-1)}, \pm\sqrt{\beta(\rho-1)}, \rho-1)$。ヤコビ行列の固有値で線形安定性が決まる。
- 分岐シナリオ: $\rho=1$ でピッチフォーク分岐(対流の発生)、$\rho \approx 24.74$ でサブクリティカル・ホップ分岐(カオスの出現)。$\rho=28$ はカオスがクリアに現れる絶妙な値。
- 初期値鋭敏性とリアプノフ指数: 最大リアプノフ指数 $\lambda_1 \approx 0.906 > 0$ がカオスの定量指標。予測可能時間は $T_{\mathrm{pred}} = \ln(\Delta/\epsilon)/\lambda_1$ で、観測精度を上げてもわずかしか伸びない。
- ストレンジアトラクタの幾何: 位相空間体積は $e^{-13.67t}$ で収縮するが、フラクタル次元 $D_{KY} \approx 2.06$ の集合に吸着される。
- Python実装:
scipy.integrate.solve_ivpのRK45で高精度に積分し、3D蝶々形を可視化。Bennettin法で最大リアプノフ指数を約 $0.90$ と推定。$\rho$ をスキャンした分岐図で理論的分岐点と数値結果の整合を確認。
ローレンツ方程式は「単純さと複雑さの境界」に位置する魅力的なモデルです。たった3本の式の中に、決定論と予測不能性、対称性の自発的破れ、ストレンジアトラクタ、フラクタル幾何、エルゴード性といったカオス理論のほぼ全てのキーコンセプトが詰まっています。一度この系を自分の手で動かして可視化してみると、気象から脳波、金融まで広がるカオス科学の世界観が、確かな実感とともに開けてくるはずです。
次のステップとして、以下の記事もご参考ください。