カオス発見の第一報を読む【論文紹介】#12
前略(あけおめ)
新年っぽいネタということで、初モノで初心に還るという意味も込めて古典的名論文の紹介をしてみたい。
Edward N Lorenz著
Deterministic nonperiodic flow
決定論的で非周期な流れ
Lorenz63.pdf (bu.edu)
”カオス(chaos)”とは論文には書いてないが、現在カオスと呼ばれる概念を最初に科学的に記述し、その美しく深遠な数学的構造を調べ、研究発展の可能性も明晰に記した名論文である。その内容の講演会を行うにあたって著者がタイトルを決めかねていたところ、主催者が『予測可能性：ブラジルの蝶の羽ばたきがテキサスに嵐を引き起こすか？(原題 Predictability: Does the Flap of a Butterfly's Wings in Brasil Set off a Tornedo in Texas?)』と勝手に決め、バタフライ効果という言葉が生まれるきっかけともなった名論文でもある(ちなみに、論文中にbutterflyという単語は一度も出てこない)。
1963年3月に発表され、60年と10カ月なのでこの際(間に合ってないとは言わない)、カオスはなんとなく知っているけど数式では全然知らないという人も視野に入れて紹介記事を書いてみたい。
ザルツマンモデルと呼ばれる次のような状況を考える(本論文では5章。4章までは今となっては前提知識のような記述で、本記事では追々説明する)。
高さHの箱の中に流体が入っている。下面と上面の温度差が一定で下面が高温とすると、下面で熱膨張した流体が浮力を得て、粘性による抵抗力を打ち破ると上昇し、対流が発生する。温度差をうまく調整すれば、ロール状の対流が発生し、y軸方向の流れはゼロで、xz平面だけで回転する流れを考えればよい。
流体一般の挙動を示すナビエストークス方程式に、ザルツマンモデルの制限を代入して整理すると次のようになる(式番号は本論文と同じ)。
$${\frac{\partial}{\partial t}\nabla^2ψ=-\frac{\partial(ψ,\nabla^2ψ)}{\partial(x,z)}+ν\nabla^4ψ+gα\frac{\partial θ}{\partial x}}$$, (17)
$${\frac{\partial}{\partial t}θ=-\frac{\partial(ψ,θ)}{\partial(x,z)}+\frac{ΔT}{H}\frac{\partial ψ}{\partial x}+κ\nabla^2θ}$$, (18)
ψ(x,z)は流線関数で、各(x,z)地点での流速を表す。θは対流がなく温度勾配が一定の場合の温度、gは重力加速度、aは熱膨張係数、vは動粘性率、κは熱伝導率を表す。ロール状の対流が定常的に起こる場合は、解を次のように置くことができる。
$${ψ=ψ_0sin(\frac{πax}{H})sin(\frac{πz}{H})}$$, (19)
$${θ=θ_0cos(\frac{πax}{H})sin(\frac{πz}{H})}$$, (20)
sin, cosによってきれいな円形の流れを表していると、ここではわかればよい(余力のある人は(19)(20)を(17)(18)に代入して成り立つと確認すること)。
ザルツマンによる先行研究で、レイリー数$${R_a=\frac{gαH^3ΔT}{νκ}}$$が次の値を超えるとロール状の対流が起こることは分かっている。
$${R_c=\frac{π^4(1+a^2)^3}{a^2}}$$
レイリー数は、熱対流が起こる力と粘性の抵抗力の比を表す無次元量で、容器の高さHや流体の種類によって粘性νが違っても、この比さえ合っていれば同じ熱対流が再現できる。Rcの最小値は、$${a^2=\frac{1}{2}}$$のときで、$${R_c=\frac{27π^4}{4}}$$を取る。レイリー数は、上下の温度差ΔT以外は実験中には変えられない定数なので、温度差がある値より大きければ対流が発生するということである。
ザルツマンは、(19)(20)の$${ψ_0, θ_0}$$を求めるためにψとθを二重フーリエ展開して、
$${\frac{aψ}{(1+a^2)κ}=X\sqrt{2}sin(\frac{πax}{H})sin(\frac{πz}{H})}$$, (23)
$${\frac{πR_aθ}{R_cΔT}=Y\sqrt{2}cos(\frac{πax}{H})sin(\frac{πz}{H})-Zsin(\frac{2πz}{H})}$$, (24)
の解を得た。ここで係数X,Y,Zは無次元化した時間(実験ごとにHやκが違っても、この倍率にすれば同じ速度で再生されて見える時間単位)$${τ=\frac{π^2(1+a^2)κt}{H^2}}$$のみに依存し、次のような関係を満たす。
$${\dot{X}=-σX+σY}$$, (25)
$${\dot{Y}=-XZ+rX-Y}$$, (26)
$${\dot{Z}=XY-bZ}$$, (27)
ここで、$${σ=\frac{ν}{κ}}$$はプラントル数、$${r=\frac{R_a}{R_c}}$$、$${b=\frac{4}{1+a^2}}$$である。
Xは、ロール状の対流の強さに比例、
Yは、ロール状の対流の上昇流と下降流の温度差に比例、
Zは、鉛直方向の温度勾配が一定の場合からの鉛直温度分布のひずみに比例する。
ここからは、X, Y, Zで表される3次元空間を考える。この空間中の1点が、ロール状の対流の1つの状態を表し、X, Y, Z空間中の軌跡が状態の変化を表す。(この、ロール状の対流が発生する実空間から要素を抽出した(物理学で言う)位相空間という考えが本論文2章の内容であり、今後は位相空間での振る舞いについての議論が続く。)
25~27式のドットがついた$${\dot{X}, \dot{Y}, \dot{Z}}$$はそれぞれ時間τでの微分で、ある(X, Y, Z)地点においてその値を右辺に代入し、その地点ではこの速度$${v=(\dot{X}, \dot{Y}, \dot{Z})}$$で位相空間中を移動するという意味である。そのような移動をして、位相空間中に軌道が描かれる。
位相空間中の軌道については、次のような重要な性質がある。
まず、状態を表す点が(25)~(27)のような速度場に従って運動する場合、軌道は決定論的で、位相空間中の1点を取れば、その点を通る軌道は一意に、前後ともいくらでも長く定まる。
また、他の軌道やすでに通った軌道と交わることはない。もし交わった場合、交わった点からは軌道が2つに分かれていることになり、その点からの軌道は一意に定まらず、前述の性質に反するからだ。
ちなみに、位相空間での振る舞いは、この動画シリーズの第2回と第4回がわかりやすい(ついでに第7回も見ておくとこの先の理解にも役立つと思う)。ChaosⅡ ベクトル場 （レゴの競争） (youtube.com)
(25)~(27)式を行列で書くと次のようになる。
$${\begin{bmatrix}\dot{x_0}\\\dot{y_0}\\\dot{z_0}\end{bmatrix}=\begin{bmatrix}-σ&σ&0\\(r-Z)&-1&-X\\Y&X&-b\end{bmatrix}\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}}$$
時間微分とは、微小時間の変化のことと考えられるので、右辺で、位相空間中のある点の座標$${(x_0, y_0, z_0)}$$にこの行列が1回かけられるごとに微小時間1ステップ進むと考えられる。
微小体積の変化は、ベクトル場の発散として次のように考えられる。
$${\frac{d\dot{X}}{dX}+\frac{d\dot{Y}}{dY}+\frac{d\dot{Z}}{dZ}=-(σ+b+1)}$$
これは位置や時間に問わずマイナスなので、位相空間中にある有限体積の領域を取ると、(25)~(27)式の時間発展によってその体積は常に小さくなっていくことがわかる。時間を十分に発展させれば、位相空間中のどの点も、体積ゼロの領域(平面か線か点かはわからないが)に押し込められると推測できる。
このベクトル場の固定点(状態がその点に置かれたら、速度は3成分ともゼロになり動かなくなる点)を考えると、まず(x, y, z)=(0, 0, 0)で、対流が全くない状態を表す点がある。
この場合の行列の特性方程式は次のようになる。
$${[λ+b][λ^2+(σ+1)λ+σ(1-r)]=0}$$, (32)
r>1のときこれの解は、原点ともう2つ、
(X, Y, Z)=$${(\pm\sqrt{b(r-1)}, \pm\sqrt{b(r-1)}, r-1)}$$である。
32式を展開すると
$${λ^3+(σ+b+1)λ^2+(r+σ)bλ+2σb(r-1)=0}$$, (33)
これは1つの実数解と2つの共役の複素数解を持つ。
$${λ^2}$$とλの係数の積が定数項と等しい場合は、2つの複素数解は純虚数になり、その条件は
$${r=\frac{σ(σ+b+3)}{σ-b-1}}$$, (34)
これが定常対流の不安定性に対する r の臨界値である。
σ<b+1 ならば rの正の値は (34) を満足せず, 定常対流は常に安定であり、原点以外の2つの固定点は、位相空間から全ての軌道が集まる体積ゼロの領域に含まれる。
σ>b+1 ならば定常対流は十分高いレイリー数に対して不安定になる。(34)の複素根の存在は、原点以外の2つの固定点から少しでもずれると、固定点からさらに離れる向きの速度を得る不安定な固定点となることを示す。
σがb+1を超えるかどうかで原点以外の2つの固定点は周囲の軌道を集めたり反発したりする。
ここまでの議論をまとめると、位相空間中のすべての点は最終的に体積ゼロの領域(平面か線か点かはわからないが)に押し込まれるのだが、その領域の近くに発散する点もあるということである。この一見すると二律背反のような状況だが、実際のところ軌道が最終的に閉じ込められる部分空間はどうなっているのだろうか？線形代数や解析的な理論ではこれ以上は明らかにはならなさそうなのでコンピュータの力を借りることにする。
さてここからは、(25)~(27)式の時間発展をコンピュータにやってもらう。
(本論文4章は計算の手法。1963年当時は"コンピュータを使うこと"すら珍しく、今や主流であるルンゲクッタ法もまだなく、前進差分法と中心差分法の平均を使っている。ただし、ルンゲクッタ法を使ったとしても以降の議論は変わらない。)
ザルツマンの先行研究に倣って、σ=10, a^2=1/2, b=8/3とする。
このとき、ロール状の対流が起こる最小のレイリー数は、r=470/19≒24.74なので、それより大きめのr=28でシミュレーションする。固定点は原点と点C: (X, Y, Z)=(6√2, 6√2, 27), 点C': (X, Y, Z)=(-6√2, -6√2, 27)の3つである。
初期条件としては、 対流のない状態から少し離れた (X, Y, Z)=(0,1,0) を選んだ。まず、対流の上昇流と下降流の温度差に比例するY座標の時間変化だけを見てみる。
1000ステップごとに区切っている
始め35ステップ目でピークに達しているのは、対流でかき混ぜられてない初期の上面と下面の流体が初めて出会うところで温度差が大きくなっているだけ。それ以降は、振動がありつつも定常的な対流がしばらく続く。1000ステップ以降は振幅がだんだん大きくなるのが目立ってきて、1650ステップあたりで符号が反転する。符号反転は、ロール状の対流の向きが右回りか左回りか反転することに対応する。それ以降、1~3回ピークを見せた後に符号反転する挙動を繰り返すが、ピークが何回で符号反転するかはこれだけではよくわからない。
Y座標だけでなくXYZ空間での時間発展を描いたのが次の図である。
Z軸を上方向とすると、上図は真横(X軸方向)からみた図、下図は真上からみた図
ステップ数は上2桁のみ表示
原点以外の2つの固定点C, C'の周りを周るループが見て取れる。固定点C, C'が2つのループの中心にあるということはザルツマンによる先行研究でわかっている。
Cの周りのループとC'の周りのループは、先ほどと同様にロール状の対流で言うと、右回りと左回りの対流に相当し、C周りとC'周りが移り変わるのは対流の向きが反転することを表している。
どちらのループも、外側から上を通って真ん中あたりに落ち着き、真ん中の下のほうで、次に右回りになるか左回りになるかの分かれ道がある。
X軸方向からみた図 面上の細い線はXの値の等高線 点線は重なっている部分の奥のほう
図3は、初期位置をいくつか変えて6000ステップまで計算を行い、軌道が最終的に収束する領域を描いたものである。今やよく知られたローレンツアトラクタが初めて描かれた図だ(本論文にローレンツアトラクタとは書いてないが)。C周りを周った軌道とC'周りを周った軌道は、外側から上を周って真ん中あたりで合流している。
同じ平面内でC周りとC’周りに行ったり来たりするということは軌道が交わっているのではないか、軌道の一意性に反するのではないか思われるかもしれないが、どの軌道も平面に漸近するだけで、限りなく近づくことはあっても交わることはない。合流するあたりであっても、位相空間中のある1点を取ればその前後の軌道が一意に定まることに変わりはない。
C周りかC'周りか、次にどちらになるかをよく調べてみると、zの極大値によって、次にCの周りのループかC'の周りのループに入るかが決まるように見える。それを調べたのが図4で、軌道上でZが極大値を取ったときの値を取り、n番目の極大値Mnとn+1番目の極大値Mn+1の関係をプロットしたものである。
C周りに行くかC'周りに行くかは気まぐれなように見えたが、2つの曲線にきれいに乗っている。25~27の支配方程式を知らなくても、ある時点でのZの極大値がわかれば、その次のZの極大値はかなり正確にわかるということだ。
385付近で鋭いピークが見えるが、これは図3のアトラクタで言うと、次に右に行くか左に行くかの分かれ道でY=0の境界に近く、ギリギリ下のほうまで行ってからアトラクタの外側ギリギリのところから上に行くので、Mnが385付近でMn+1が極端に大きくなる。
ここで、図4の振る舞いを少し単純化したテント写像を考える。
数列のn番目(横軸)とn+1番目(縦軸)の関係
テント写像は数式で書くと次のようになる。
$${M_n<1/2}$$のとき$${M_{n+1}=2M_n}$$
$${M_n=1/2}$$のとき$${M_{n+1}=1}$$
$${M_n>1/2}$$のとき$${M_{n+1}=2-2M_n}$$
テント写像に従う2つの数列$${M_0, M_1, M_2…}$$と$${M'_0, M'_1, M'_2…}$$を考える。
初期状態はわずかな差εしかなく、$${M'_0=M_0+ε}$$とする。
n番目の項は$${M’_n=M_n\pm2^nε}$$と表される。
テント写像の表式を見ればわかる通り、誤差εは、1回写像するたび2倍になる。Mnが1/2を超えるかで、Mnにかけられる係数は2か-2で符号は違うが、誤差は$${\pm2^nε}$$として指数関数的に開いていく。
これが、カオスの本質である初期値鋭敏性である(本論文では「軌道の不安定性(unstability)」と書いているが)。
図4でも同様に、ある時点でのZの極大値を測定しても、もしそこにわずかな測定誤差があれば、その差は指数関数的に開いていく。このように指数関数的に差が開くことは本質的に流体の現象として備わっていることである。
図4は図5のテント写像と比べてやや曲がっているが、その曲がりは不安定性には関係ない(このことが数学的に厳密に証明されたのは2001年のことだが)。
図2の位相空間で考えても、すぐ近くの軌道、それこそ2つのロール状対流の状態の差が蝶の羽ばたき程度の風速の差でも、2つの軌道の差が指数関数的に開いていき、それなりに時間が経てばトルネード1つ分くらいの大きな差になることは十分にあり得る。
ここまでで、ザルツマンが考えたロール状の対流に見られるカオス的振る舞いを見てきたが、他の流体系にも同様のことが言えるだろうか？
気象現象の長い観測史の中で、大気がぴったり周期的であることは今まで観測されていない。類似現象が起きたとしても、その後しばらくすれば全く異なった状態になる。これはつまり、どのような予測方式もある程度の短期間なら予測できるが、長期間では”必ず”正しい結果を与えることはできないことになる。
ではその長距離とはどのくらいなのかという問題は、本論文の結果からはわからない。「おそらくそれは数日かもしれないし、数世紀かもしれない」と本論文には書いてあるが、現在の日本での気象予報技術では、10日くらいが限度とされている。また、これは気象庁しか公開していないが、週間予報には信頼度があり、数日先でも信頼度が低かった部分は予報が変わったりする。
本論文でのロール状対流のような理想化された系であれ、大気をできるだけ忠実に再現するシミュレーションであれ、初期条件を少しずつ変えた数値解の組を比較することによって答えを得ることができる。それによって、本論文ではバタフライ効果が見出され、現実の気象予報の場合、軌道のばらけ具合が(おおむね)予報信頼度となる。(本論文ではこうは書いてないが、ほぼこれに相当することが書いてあって、1963年とは思えない先見の明に驚く。)
以下、私の感想だが、改めて見るとこの論文の凄さは、コンピュータの力に頼ってすごい結果を得ただけでなく、考察の鋭さとその影響の大きさだと思う。
コンピュータは、計算をぶん回すだけでは、解析的に出せない複雑な方程式の解を出す道具でとどまっていただろう。そこに、Lorenzの鋭い考察によって、カオスという現象が見出され、しかもそれが気象をはじめとして身の回りにあふれていて様々な分野で同様の現象が発見され、一つの新しい科学の分野が作られたと言ってもいい。そして今や、カオス現象を探索するためにコンピュータが投入されることも多々ある。
近年、コンピュータのパワーは圧倒的になり、LLMや汎用人工知能の発展によって、人間の知能は相対的に矮小なものになっていくだろう。幹線道路をびゅんびゅん車が走っている横をとぼとぼ歩いているような存在になるだろう。しかし、幹線道路からでもちょっと脇道に逸れれば車ではたどり着けない見たことのない景色が開けているかもしれないし、新しい道を敷く工事をして新たな知の地平を切り開いていけるかもしれない。人間の知能はそのようなものでありたいとこの論文を読んで思った。
この記事が気に入ったらサポートをしてみませんか？