Analysis of Heteroclinic Turbulence in Nonlinear Reaction Diﬀusion Equation

(1)

非線形反応拡散方程式におけるヘテロクリニック乱流の解析

Analysis of Heteroclinic Turbulence in Nonlinear Reaction Diﬀusion Equation

早稲田大学大学院理工学研究科物理学及応用物理学専攻

統計物理学研究

渡橋憲司

2010 年 6 月

(2)

序文

乱流とは時間・空間的な乱れが持続する現象を指し、組織だった乱流研究としては 1883年にＲｅｙｎｏｌｄｓの行った、円管を内臓した水槽を用いた流れの観測実験を発端とする。その後、コンピューターの誕生・発達に伴い、1970年代からナヴィエ・ストークス方程式のシミュレーションを用いた乱流研究が目覚ましい発展を遂げた。流体乱流の解析とともに、シュミレーションを用いて、流体以外の系を記述する様々な非線形偏微分方程式にも乱流が存在することが明らかになり、乱流研究の対象は広がっている。

シミュレーションにより得られた乱流解に対する解析では、主に統計手法が行われてきた。具体的にはスペクトル解析や相関解析が挙げられる。Ｋｏｌｍｏｇｏｒｏｖの、等方性乱流の空間スペクトルの5/3乗則を導く現象論的解析に始まり、スペクトルが数値計算により具体的に求められるようになると、様々な方程式の乱流解について調べられた。一方、レイノルズ数の研究の流れから、系の特性長・特性時間についても研究がなされ、またカオス理論の発展に伴い乱流のフラクタル性やリャプノフ指数についての数値的研究も有力な手法となっている。そうした統計解析の発展研究のひとつに、統計力学と臨界現象の概念を用いた解析がある。臨界現象の理論を乱流のフィールドに適用し、系の特性量として相関長・時間やリャプノフ指数を選び、それらについて様々なパラメータ依存性を調べ、乱流の普遍性についての知見を得ようとする試みがＨｏｈｅｎｂｅｒｇ等によりなされ、乱流の理解を大きく進展させた。

反応拡散方程式の乱流解の発生点付近の振る舞いについての有効な理論的解析として、位相縮約法が挙げられる。反応拡散方程式の解の分岐点近傍での縮約方程式を用いて乱れの発生時に関する理論的な知見を得ることが可能となったが、乱流の大域的構造については今なお未解明の部分が多い。

最近見つかった乱流の新たな発生機構として、生態系（ロトカ・ヴォルテラ系）や遺伝子発現系（レプリケータ系）に典型的にみられる、ヘテロクリニックサイクルによる乱流現象がある。ヘテロクリニックサイクルとは、複数の不安定な鞍部点とそれら

(3)

をつないだ軌道のことで、そのサイクル付近の解は、サイクル自身ににたどり着かずに巡回しつづけ、その鞍部点付近の滞在時間は指数関数的に長くなる特徴を持つ。ヘテロクリニックサイクル近傍の解が拡散項の効果により乱流化することが2003年にペトロヴスキーによって明らかにされたが（この現象をヘテロクリニック乱流と名付ける）、その乱流解の振る舞いについては、ほとんど明らかになっていない。

そこで本論文では、ヘテロクリニック乱流という未知な乱流現象に対し、具体的な反応拡散方程式の乱流解のシミュレーションを行い、その発生メカニズム・特性を明らかにすることを目的とする。

(4)

1 序論

1.1 乱流・時空カオス

乱流とは時間・空間的な乱れが持続する解を指し、自然界では昔から水の流れの中で観測されていたが、最初に行われた組織だった研究は、実験的には1883年のReynolds の行った、円管を内臓した水槽を用いた流れの観測実験であり[1]、理論解析としては

1925年にPrandtlの行った、気体の分子運動論との類推で平均自由行路に類似する混合

長という概念を用いた混合長理論であった[2]。その後理論の分野では、Taylorが1935 年にはじめた等方性乱流の取り扱いに統計的手法を用いた乱流の統計理論を中心に進

んでいく(例：Kolmogorovの等方性乱流に対する現象論的理論による、空間スペクト

ルの5/3乗則)[3, 4, 5, 6]。

実験・理論と異なる第3の研究として、のちに乱流研究の主流にまで発展する偏微分方程式にシミュレーションを用いた解析もこの時期に誕生した。偏微分方程式のシミュレーション解析は、1910年のRichardsonの、楕円型方程式に対する反復法の研究、に端を発する[7]。この研究では、当時まだコンピュータが実在してなかった為、

人を雇い直接計算させた結果を発表している。その後、様々な差分法の提案（FTCS

（Forward Time Center Space）法、リープマン法[8]、蛙跳び法[9]、クランク・ニコル

ソン法[10]、2段階ラクス・ヴェンドルフ法[11]、アダムスバスフォース法[12]など）

や差分法の安定性（ノイマンの安定性解析、クーラン・フリードリヒ・レヴィーの安定性解析（ハートの安定性解析）[13, 14]など）についての研究が進み、数値計算を行う上での土台が築かれた。コンピュータの発明・発達に伴い、1960年代後半からナヴィエ・ストークス方程式の乱流解のシミュレーションによる研究が本格的に行われるようになり、乱流研究が飛躍的に進む[15, 16, 17, 18, 19]。

乱流研究と独立にかつ時を同じくして、カオス現象が発見された。カオスとは、決定論的な系に現れる一見無秩序だがその背後にいくつかの秩序構造を持つ振る舞いを指し、初期値鋭敏性を本質的に持つ力学的現象である。また非可積分であるため、シ

(8)

ミュレーション解析が主に用いられる。1899年にPoincareが3体問題では解析解が得られないことを証明し[20]、カオスの存在を示唆していたが、具体例としては、1961 年に上田による電気回路の引き込み現象を表す非線形常微分方程式の不規則現象の研究や、1963年にLorenzによるローレンツモデルに現れる非周期的な流れの研究[21]、

そして1976年にMayの行ったロジスティック写像の研究[22]で発見された。カオスについては多くの研究がなされ、対象は多岐に渡り[23]、その中で乱流を時空カオスと呼び、カオス現象の一つとしてとらえる見方が生まれ、乱流研究の大きな分野の一つとなった[24, 25, 26, 27]。

1.2 様々なモデルの乱流

コンピュータの発達によるシミュレーション研究の進歩により、乱流の主な研究対象であった流体系以外でも、様々な偏微分方程式について乱流解が存在することが明らかになってきた[28, 29]。複数の非線形偏微分方程式についての乱流解のシミュレーション結果を、乱流の発生機構に焦点を絞って挙げる。

1.2.1 ナヴィエ・ストークス方程式

粘性応力が流れに垂直な方向の速度こう配に比例するというニュートンの粘性法則を運動方程式に導入する取り組みが1882年にNavierによってはじめられ[30]、1885年

Stokesによって大成された[31]。これが今日ナヴィエ・ストークス方程式と呼ばれる粘

性流体力学の基礎運動方程式である。

ナヴィエ・ストークス方程式は、非圧縮流体の速度場の時間発展を記述する方程式で、下の形で表される。

d_tu=∂_tu+ (u∇)u=−1

ρ∇p+ν∇²u+f (1.1) ここで、u≡u(x, t)は連続体の速度場で、p≡p(x, t)は圧力、f ≡f(x, t)は外力、ρは密度（非圧縮なので一定値）、そしてνは動粘性係数である。この式は、流体要素が圧

(9)

力、粘性力、そして外力によって加速されていることを表す。この方程式で右辺の第 2項（粘性力）を省略した方程式をオイラー方程式と呼ぶ。

(1.1)式の流れの代表的な長さスケールをL、代表的な速さスケールをUとすると、

これらを用いてレイノルズ数RはR= (U L)/νと定義される。この方程式はレイノルズ数が低い時は層流だが、高い時は微小な乱れが(u∇)uの移流項の非線形性により拡大し、乱流化する。また、その乱流化のメカニズムは非常に複雑な過程をたどる。例えば、直方体容器の中のベナール対流の系では、容器の形やプラントル数（動粘性係数と温度拡散係数の比）の差によって、レイノルズ数を上げると定常解→周期運動→2 重周期運動→· · ·→n重周期運動→· · ·→乱流というように遷移するランダウ・ホップのシナリオや、定常解→周期運動→2重周期運動→3重周期運動→乱流というように遷移するリュエル・ターケンスのシナリオ、定常解→乱流というようにいきなり遷移するシナリオ（ローレンツ系に似た遷移）、そしてランダウ・ホップのシナリオに類似しているが、周期が次々と半分になる低周波分岐、など様々な分岐過程が起こる[32]。

また、定常外力が駆動するナヴィエ・ストークス方程式では、レイノルズ数を大きくするに従って、定常解→周期解→2重周期解→3重周期解→2重周期解→乱流→周期解→2重周期解→乱流→· · ·という、非常に複雑な遷移過程をたどる場合もある[33]。

1.2.2 複素ギンズブルグ-ランダウ(cGL)方程式

複素ギンズブルグ-ランダウ(cGL)方程式は物理学者の間で最も研究された方程式の一つで、以下の形をしている場合もある。

∂_tA=A+ (1 +ic₁)∇²A−(1−ic₃)|A|²A (1.2) Aは時間tと空間⃗xの複素関数で、理想的な自励振動場を表し、c1 とc₃ はそれぞれ線形と非線形散逸のパラメータである。cGL方程式(1.2)は現象論的には Newell と Whiteheadによって1971年に[34]、また化学系の分野でKuramotoとTsuzukiによって1976年に導出された[35]。この方程式は反応拡散方程式のホップ分岐近傍での縮約方程式として導出することが出来る。

(10)

トラヴェリングウェーブA_k(x) =a_kexp^i(kx⁻^ωt)がすべての波数について不安定になる条件、つまり1−c₁c₃ <0の際に系が不安定になり（=ベンジャミン・フェアー不安定）、乱流が生じることが知られている[36]。乱流はベンジャミン・フェアー不安定のライン付近では、位相乱流と呼ばれる、振幅|A|はほぼ1がその位相要素が乱流的に振る舞う弱い乱れが現れる。その際の特徴として、ワインディングナンバーが初期状態を保存し続けることと、相関が非常に長いことが挙げられる。但し、ワインディングナンバーが位相乱流領域で本当に初期状態を保持しつづけるかどうかは未だ議論の余地を残している[37]。

一方、ベンジャミン・フェアー不安定ラインから遠く離れたパラメータでは、振幅

|A|はもはや、1から大きく揺らぐ、振幅乱流と呼ばれる強い乱れが生じる。振幅乱流では欠陥（ディフェクト）と呼ばれる、局所的に振幅が0になる点（位相部分に不連続が生じる点）が現れる。他にも、パラメータによっては、振幅乱流と位相乱流の間に、バイカオスと呼ばれる、両方の乱流が共存する領域が生じたり、間欠的な乱流の発生する領域が出現する[38]。

空間2次元の場合は、振る舞いはさらに多様になる。点欠陥（ポイントディフェクト）に関して、位相乱流や振幅乱流が現れるだけでなく、標的パターンや回転らせん波のパターンが現れることが知られている[39]。

1.2.3 蔵本・シバシンスキー(KS)方程式

蔵本・シバシンスキー(KS)方程式は、cGL方程式の位相乱流領域で、位相成分ψを近似することにより、以下のような形で得られる。

∂_tψ =−∇²ψ− ∇⁴ψ+ (∇ψ)² (1.3) またはv = 2∇ψを用いて、

∂_tv =−∇²v− ∇⁴v+v∇v (1.4)

(11)

という形になる。この方程式は、1984年にKuramotoによって提案された[40]。またそれより以前に層流の炎の表面の振る舞いを記述する方程式として1977年にSivashinsky によって提案された[41]。この方程式のパラメータは、システムサイズのみである。システムサイズを増やすと、不動点から、単峰型の不動点、周期軌道、2つの山を持った不動点、振動解と乱流のマルチベイシン、3つの山を持った不動点、振動解と乱流のマルチベイシン、4つの山を持った不動点、と遷移過程を経て、乱流のシングルベイシン領域へと移り変わる[42]。その際、周期倍加分岐が起こる[43]。

1.2.4 グレイ・スコットモデル

グレイ・スコットモデルは3次の解放系での自己触媒モデルとして1984年にGray

とScottによって提案され[44]、以下のような形をとる。

∂_tu=D_u∇²u−uv²+F(1−u)

∂_tv =D_v∇²v+uv²−(F +k)v

(1.5)

ここで、uは反応物Uの濃度、vは自己触媒V の濃度、F は外部からのU の流入率、

F +kはV の流出率である。また、自己触媒反応はU + 2V →3V、V →Wという式であらわされる。この系では、乱流がホップ分岐を経て発生し、その後、安定な不動点と不安定な不動点、そして空間非一様な定常解と、それぞれの解に漸近するが、どの解にも落ち切らず、巡回を続ける、ヘテロクリニックサイクルによる乱流が発生する[45, 46, 47]。

1.2.5 フィッツヒュー・南雲モデル

フィッツヒュー・南雲モデルは、神経細胞の興奮・抑制状態や、化学反応を記述する方程式として用いられる。また、白金上での一酸化炭素の酸化現象を記述する方程式でもあり、以下の形を持つ。

∂_tu = −1

ϵu(u−1)(u−b+u

a ) +∇²u

(12)

∂tv = f(u)−v (1.6)

f(u) =











0, 0≤u <1/3,

1−6.75u(u−1)², 1/3≤u≤1, 1, 1< u,

ここでuは一酸化炭素の白金上での吸着量、vは表面での割合、aは酸素の圧力、bは一酸化炭素の圧力、そしてϵは温度を表すパラメータである。この系でも、系は乱流解をもつ。その際の乱流化のメカニズムは、ホップ分岐により乱流になり、その後ヘテロクリニックサイクルによる乱流が発生し、安定な不動点と不安定な不動点を巡回しつづけるヘテロクリニックサイクルがトラヴェリングウェーブ化するまで続く[48]。

1.2.6 ニコラエヴスキーモデル

ニコラエヴスキーモデルは、粘弾性媒質での非線形効果も考慮にいれた、振動の時間発展を記述する方程式で、Nikolaevskiiによって1989年に提案され[49]、以下の形を持つ。

∂_t=−∇²[ϵ−(1 +∇²)²]ψ−(∇ψ)² (1.7) ここでψは媒質の振動、ϵは分岐パラメータである。よってこの系のパラメータは、L とϵの2つである。この系は振動場のチューリング不安定による乱れとして、縮約法により導出される[50]。ϵの変化に対し、ϵ = 0でスーパークリティカル分岐を起こし、

ϵ >0で乱流が発生する[51]。またLの変化に対し、L >2π/(1 +ϵ^0.5)^0.5で乱流が発生するが、その分岐はT²トーラス倍加分岐である[43]。

1.2.7 バーガース方程式

バーガース方程式は、振動理論や、音響学、プラズマ物理などで様々な非線形波を記述する方程式として知られており、下のような形を持つ、

∂u

∂t =−ν∂u

∂x +v∂²u

∂x² (1.8)

(13)

ここで、uは非線形波の速度場を表す。この方程式はBurgersによって乱流のモデルとして1939年に提案された[52]。他の散逸系の方程式とは異なり、可積分系であり、ホップ・コール変換により解けることが知られている[53, 54]。具体的には

u=−2ν 1 U

∂U

∂x (1.9)

の変換を用い新たな変数U による表示を使うと、バーガース方程式は

∂U

∂t =ν∂²U

∂x² (1.10)

という拡散方程式に帰着出来、その結果、解は u(x, t) = −2ν ∂

∂xln[(4πνt)⁻^0.5

∫ _∞

−∞

exp{−(x−y)² 4νt − 1

2ν

∫ _y

0

u(z,0)dz}dy] (1.11) と表される。その為、カオスの性質（初期値鋭敏性）を持たない。具体的には初期値として乱れた形を代入すると、有限時間、乱れが減衰しながら振動を続け、無限遠で完全に減衰する、という振る舞いを見せる。

この方程式は可積分系なので、相関やスペクトルといった乱流によく用いられる統計について理論的な計算が可能である（詳細は本章のスペクトルの節で触れる）。

1.3 _{乱流の解析手法}

乱流を解析する確立された理論はないように思える。ただこれらの状況下でも、乱流を解析する手法がいくつか提案された（すぐれたまとめとしてCrossとHohenberg(1990)[55]

を見よ）。

1.3.1 スペクトル解析

空間的広がりを持った平均値0の時系列a(x, t)について、ある時刻tでの空間スペクトルS(k;t)は次のように定義される。

S(k;t) =| 1

√2π

∫ _∞

−∞

a(x, t) exp (−ikx)dx|² (1.12)

(14)

またある位置xでの時間スペクトルS(ω;x)も同様に、

S(ω;x) =| 1

√2π

∫ _∞

−∞

a(x, t) exp (−iωt)dt|² (1.13) と定義される。時系列の平均値が0でない場合は、a^′ =a− < a >と時系列を変数変換すれば良い。乱流でこれらの量が求められる際は、空間スペクトルについては時間t の平均を施したS(k)が、時間スペクトルについては空間xの平均を施したS(ω)について、計算されるのが通常である。

スペクトルは、フーリエ級数・変換を用いた研究から生まれた。フーリエ級数は、

1811年にFourierが熱伝導の問題の解法に用いたのが最初である[56]。

空間スペクトルについて、いくつかの方程式についての結果を述べる。ナヴィエ・

ストークス方程式で等方性乱流については、空間3次元では、低波数領域ではエネルギーピークを持ち、中波数領域で-5/3乗のべきを持つ慣性領域を持ち、高波数領域で指数減衰する[5]。一方空間2次元では、低波数領域で-5/3乗のべきを持ち、高波数領域では-3乗のべきを持つことが分かっている[57]。KS方程式では、空間1次元空間スペクトルは低波数領域ではフラット、中波数領域でピークを持ち、高波数領域で指数

減衰する[58]。また、ゼロレベルクロシングの空間スペクトルはべき減衰する[59]。複

素・ギンズブルグランダウ方程式については、位相乱流については、KS方程式と同じスペクトル形を持ち[38]、振幅乱流については、低波長領域でフラット、中・高波数領域で指数減衰する[60]。ニコラエヴスキーモデルでは、分岐点近傍(ϵが小さい正の値) では、低波数領域でフラット、中波数領域で一度減衰後ピークを持ち、高波数領域で振動しながら0へ減衰する値を持ち、分岐点から離れた場合(ϵが十分大きい値)、KS 方程式と同様のスペクトル形を持つ[61]。バーガーズ方程式については、可積分系であるため、初期状態に依存して空間スペクトルが決まる。また、限られた場合では理論的に計算することが可能である。初期速度場が一様にランダム挙動を行う、つまり、

スペクトルが

S₀(k) = α²|k|ⁿb₀(k) (1.14) で表されると仮定する。ここでnは初期状態を決定するスペクトル指数で、b0(k)は、

(15)

b0(0) = 1でb0(k)→0 (k → ∞)に連続的に減衰し、その減衰はどんなべきよりも速いと仮定する。その時、空間スペクトルは、

S(k)∝











|k|⁻² (|k|>> k₁の時)

|k|² (k₂ <<|k|<< k₁の時)

|k|ⁿ (|k|<< k₂の時)

(1.15)

となることが知られている[62, 63]。ここでk₁、k2はnと時間によって変化する値である。

1.3.2 自己相関関数

空間的広がりを持った平均値0の時系列a(x, t)について、ある時刻tでの空間自己相関関数C(x;t)は次のように定義される。

C(x;t) = lim

X→∞

1 2X

∫ X

−X

a(x^′, t)a(x^′+x, t)dx^′ (1.16) またある位置xでの時間自己相関関数C(t;x)も同様に、

C(t;x) = lim

T→∞

1 2T

∫ T

−T

a(x, t^′)a(x, t^′+t)dt^′ (1.17) と定義される。乱流でこれらの量が求められる際は、空間自己相関関数については時間tの平均を施したC(x)が、時間自己相関関数については空間xの平均を施したC(t) について、調べられるのが通常である。

ウィナー・ヒンチンの定理より、空間自己相関関数は空間スペクトルのフーリエ変換、時間自己相関関数は時間スペクトルのフーリエ変換として、それぞれ求められることが知られている[64]。

空間相関関数についての結果は、1次元のcGL方程式の位相乱流では指数減衰し[65]、

振幅乱流では指数減衰が崩れる[66, 38]。また1次元のKS方程式では指数減衰とべき減衰が現れる[67]。そして1次元の拡散ロトカヴォルテラ方程式系などで[68]、調べられた。時間相関関数については1次元のKS方程式[69]で調べられた。

(16)

1.3.3 臨界現象の手法を用いた解析

乱流の統計解析の発展研究として有力なもののひとつに、乱流に統計力学と臨界現象の概念を用いた解析がある。臨界現象の理論では、臨界パラメータ(系の分岐パラメータ)に関連した物理量の臨界指数が重要となる。多くの指数が実験的、および数値計算的に調べられ、それら2つの結果が比較された。さらには、スケーリング仮説

(ウィドム [70]、ルーシュブルーク [71]、フィッシャー [72]、そしてジョセフスン [73]

など)、斉次性仮説そして普遍性仮説が提案され[116]、最後に、それらの仮説はウィルソンによるくりこみ群を用いた解析により、理論的に証明された [75]。これらの概念を乱流のフィールドに適用し、系の特性量について様々なパラメータ依存性を調べる手法が提案された [38]。この手法の目的は分岐点付近の臨界指数を調べ、乱流の普遍性に関する情報を手に入れることである。

その際、乱流の特性量としてなにをとるか、という問題がある。これまでの解析で、

統計学的に、もしくは力学的に以下で挙げる量が主に調べられてきた。

1.3.4 相関長

時空カオスの現れる偏微分方程式系で、空間自己相関関数(1.16) を用いて相関長ξ について、もっとも単純なものとして、次のような定義が提案された[38]。

C(r)∼exp(−x

ξ) [r→ ∞] (1.18)

この定義は指数長とも呼ばれ、1次元のcGL方程式[38, 76]、1次元のKS方程式[65]、

そして1次元の反応拡散系[68]などで調べられた。しかしこの定義には、系の相関関数が指数的に減衰することを仮定しなければならないという問題点がある。

ほかには、系の従属変数u(x)について、次のような量 K(l) = <(u−< u >_l)⁴ >_l

3<(u−< u >_l)²)>²_l −1 (1.19) を定義する[77]。ここで< A >_l=∫

Ω(l)dxA(x)である。lが相関長より充分大きい時は

u(x)が正規分布に従うので、K(l)は0となる。しかしlが相関長より小さい時はu(x)

(17)

は正規分布に従わず、K(l)は有限の値を持つ。この性質を用い、K(l)が0から有限の値に変化するlを相関長とする方法が提案された。この定義には、問題点として、数値計算的にも相関長を求めるのに非常に高いマシンの性能を必要とするということが挙げられる。

他には、第iリャプノフ指数のリャプノフベクトルu¯_i(x, t)を用いて ξ_i^d= [∫

d^dx¯u²_i(x, t)]²

∫ d^du¯⁴_i(x, t) (1.20)

というように相関長を定義する方法も提案された[78]。ξ_i^dはtの関数であるが、積分の区間を充分広くとることにより、空間平均の値もほぼ一定値になる。

1.3.5 相関時間

相関時間についても、相関長(1.18)と同様に、時間自己相関関数(1.17)を用いて C(t)∼exp(−t

τ) [t → ∞] (1.21)

のように定義される。相関時間については、cGL方程式[38]やKS方程式[66]で調べられた。

1.3.6 リャプノフ指数

乱流をカオス現象として見る際のカオスの尺度として考えられているのが、リャプノフ指数である。

第iリアプノフ指数λiは接空間上の2つの近い軌道同士の指数関数的離散率として以下のように定義される。

λi = lim

t→∞t⁻¹ln|DF|i, (1.22) ここで(DF)_ij ≡∂F_i/u_jは軌道から得られたヤコビアン行列( ˙u = F(u))、|A|iは正方行列Aのi番目の固有値、そして|A|i−1 ≥ |A|i ≥ |A|i+1とする。λiはグラムシュミットの直交化法で得られる[79]。

(18)

第1リアプノフ指数が正の場合はカオス、0の場合は周期解、負の場合は定常解を表す。また、リアプノフ指数間の特徴として、λ1 ≥ λ₂ ≥ · · · ≥ λ_nとなることが知られている。

リャプノフ指数から得られる量として、乱流の次元D_λはリャプノフスペクトルλ_i から求められる[55]。

D_λ =k−(

∑k i=1

λ_i)(λ_k+1)⁻¹, (1.23) ここでkは∑_k

i=1λ_i ≥0を満たす最大整数である。

cGL方程式についてはリャプノフ指数について調べられ、位相乱流ではほぼゼロで

あり[80]、リャプノフ次元がシステムサイズに線形に増大することや[65, 76, 81]、シ

ステムサイズが小さい時、比例関係の破れが生じることが知られている[82]。そしてリャプノフスペクトルの概形についても詳しく調べられ、リャプノフ指数が大きい時、

4つの隣同士のリャプノフ指数が並ぶことが報告されている[83]。またKS方程式についても、リャプノフ次元がシステムサイズに線形に増大し[84]、リャプノフスペクトルについても調べられた[85]。またニコラエヴスキーモデルに対して、リャプノフスペクトラの概形が調べられ、またリャプノフ次元がシステムサイズに比例することが知られている[86]。

1.3.7 縮約法

縮約法とは、系を記述する方程式の解が、アトラクターに漸近後緩慢な運動を行うような場合、緩慢な運動を記述する簡潔な方程式を導出する手法である。この方法は、

系が分岐を起こして乱雑な挙動に変化する際の、分岐点近傍で威力を発揮する。ホップ分岐の場合の反応拡散方程式について縮約法を行うと次のような結果になる。以下のような反応拡散方程式

∂X

∂t =F(X;µ) +D∇²X (1.24)

(19)

を考える。ここでXは扱う系の従属変数、Dは非負の要素を持つ対角行列、µは力学系の分岐パラメータで、力学系X˙ =Fの定常解がホップ分岐、つまり線形化方程式の固有値について、µ= 0で1対の複素共役な固有値が虚軸を横切り、µ >0で系が不安定化すると仮定する。また空間は無限の広がりをもつものとする。

このような方程式について、縮約法を分岐点近傍で実行すると縮約方程式を導出することが出来る。最終的に(1.24)式からcGL方程式

∂W

∂t =µλ₁W −g|W|²W +d△W (1.25) が導出される。各パラメータ、変数の元の系との対応関係、そして道中の計算過程については、とても煩雑な為APPENDIX Cに載せる。(1.25)式は、適切なパラメータ変換によって(1.2)式の形に直すことが出来る。ただし、(1.25)のパラメータの値によっては、直すことが出来ない場合も存在する。それについてもAPPENDIX Cで言及する。この手法は1976年にKuramotoとTsuzukiによって提案された[35]。また、詳細な解説が[40]にまとめられている。

1.4 ヘテロクリニックサイクル

ヘテロクリニックサイクルとは、2つ以上の不動点とそれらをつなげる軌道を含んだサイクルを指す。ヘテロクリニックサイクルの振る舞いは、1975年にMayとLeonard によって、ロトカヴォルテラ方程式(メイ・レオナルドタイプ) について数値計算的に発見された[87]。それとは独立に1980年にBusseが回転層内の伝送の系から、May達と同じ形の系を導出し、ヘテロクリニックサイクルの振る舞いを調べた研究により、ヘテロクリニックサイクルの研究が活発に行われることになる[88]。ヘテロクリニックサイクルの定義は、Guckenheimer達によってなされた[89]。また漸近安定なヘテロクリニックサイクル付近の解は、不動点付近に滞在しながら巡回を繰り返し、その滞在時間を指数関数的に伸ばす、という性質を持つ[88]。

漸近安定なヘテロクリニックサイクルを持つ力学系に拡散項を付けると、乱流が起こる。このタイプの乱流を今回我々はヘテロクリニック乱流と命名することにする。

(20)

ヘテロクリニック乱流の例としてグレイ・スコットモデルや、フィッツヒュー・南雲方程式などが挙げられる。グレイ・スコットモデルは偏微分方程式であるが、不安定なホップ分岐点から安定なサドルノード分岐点へ、そして空間非一様な定常パターンから元の不安定なホップ分岐点へと、巡回を続けるヘテロクリニックサイクルが存在

する[46, 45]。一方フィッツヒュー・南雲方程式では、安定な不動点と不安定な鞍部点

によるヘテロクリニックサイクルが存在する[48]。そして、滞在時間を指数関数的に伸ばすタイプのヘテロクリニックサイクルを持つ典型的な生態モデルにレプリケータ方程式とロトカヴォルテラ方程式がある。

1.4.1 レプリケータ方程式

レプリケータ方程式は個体群における戦略頻度の時間発展を記述する常微分方程式で、進化ゲームで良く知られている。個体群が頻度X₁からX_nを持つn個のタイプE₁ からE_nに分けられ、、個体群が非常に大きく、また世代が互いに連続的に重なっていると仮定すると、方程式は下のように表される。









 dX_i

dt =X_i(

∑n j=1

g_ijX_j −

∑n j=1

∑n k=1

g_jkX_jX_k)

∑n i=1

X_i = 1

ここでg_ijは戦略jに対する戦略iの得る利得を表す行列である。右辺第一項内の

∑n j=1

g_ijX_j はE_iの適応度を表し、

∑n j=1

∑n k=1

g_jkX_jX_kは平均適応度を表す。レプリケータ力学系は最初にTaylorとJonkerによって提案され [90]、SchusterとSigmundによって名づけら

れた[91]。この方程式はn人ゲームへの一般化や[92]、無限個の戦略などの拡張が提案

された[93]。3種の場合の相図は研究され[94]、また3変数以上のレプリケータ方程式

はヘテロクリニックサイクルを持つ[94, 95]。進化ゲーム理論で現れる重要な概念であるナッシュ均衡や進化的に安定な戦略の議論において、この方程式は頻繁に用いられる。すぐれたまとめとしては[96, 97]がある。

(21)

この方程式の3種の場合のヘテロクリニックサイクルは、3つの不安定な鞍部点により形成されており、漸近安定である。そのヘテロクリニックサイクル付近の解の鞍部点周りの滞在時間は指数関数的に増大する[95]。

1.4.2 ロトカヴォルテラ方程式

ロトカヴォルテラ方程式は、n種の個体数の時間発展を記述する方程式で、数理生態学の分野で頻繁に用いられる方程式で、次のように記述される。

dYi

dt =Y_i(r_i+

∑n j=1

a_ijx_j),; i= 1,2, ..., n (1.26) ここでY_iは個体群密度、riは内的増加率、aijはi個体群に対するj個体群の影響を表し、増進効果ならば正、阻害効果ならば負となる。

レプリケータ方程式とロトカヴォルテラ方程式については、次のような関係がある。

r_i =g_in−g_nn、aij =g_ij−g_nj、yi =x_i/x_n と変換すると、n−1種のロトカヴォルテラ方程式はn種のレプリケータ方程式と一致する[98]。ただし、これは常微分系でのみ成り立ち、偏微分系では成り立たない。また、3変数のロトカヴォルテラ方程式（メイ・

レオナルドタイプ）は、3つの鞍部点によるヘテロクリニックサイクルを持ち[87]、その鞍部点付近の滞在時間も指数的に伸びる[88]。

1.5 _{本論文の目的}

本章の概説の中で説明してきたように、偏微分方程式の乱流解をシミュレーションで計算し、それに対し統計解析を行うことは、乱流を理解する上で有効な手段である。

また、ヘテロクリニックサイクルを拡散項により連結する際現れる乱流の発生メカニズム・性質は、まったく未知である。さらに、滞在時間を指数関数的に伸ばすタイプのヘテロクリニックサイクルを用いた乱流には、解のヘテロクリニックサイクル付近での滞在の影響による乱れの阻害、といった他の乱流にはない特徴が現れる事が期待

(22)

される。そして、乱流の振る舞いを調べる際、複数の同様な性質を保有する系を調べ比較する事は、その性質をより明らかにする、有用な方法である。

そこで本論文では、ヘテロクリニック乱流という未知な乱流現象に対し、具体的な複数の反応拡散方程式の乱流解のシミュレーションを行い、その発生メカニズム・特性を明らかにすることを目的とする。

2章では、ヘテロクリニック乱流を実際にシミュレーションし、乱流の性質を明らかにする。具体的には空間1次元3種の拡散レプリケータ方程式を周期境界条件で扱う。

3章では、2章と異なる系についてヘテロクリニック乱流を調べ、両者の比較により、

ヘテロクリニック乱流の性質を明らかにする。具体的には競争拡散ロトカ・ヴォルテラ方程式（メイ・レオナルドタイプ）の周期境界条件下での計算結果を述べる。4章では、空間2次元のヘテロクリニックサイクルによる乱流パターンの振る舞いを明らかにする。空間2次元の3変数競争ロトカ・ヴォルテラ反応拡散方程式を周期境界条件で扱う。５章では結論と議論、そして今後の展望について述べる。

(23)

2 拡散レプリケータ方程式

本章では、ヘテロクリニック乱流を持つ系として、空間1次元で3種の拡散レプリケータ方程式を扱う。最初に、あるパラメータでの系の分岐図を、線形安定性解析と数値計算を用いて描く。それにより乱流が空間一様解のスーパークリティカルホップ分岐とともに発生することを示す。第2に、スーパークリティカルホップ分岐点付近の乱流の統計的性質を、数値計算を用いて詳細に調べる。具体的には、相関長、相関時間、第1リアプノフ指数、そしてリアプノフ次元を採用する。さらにはそれらの特性量がスケーリング則を持つことを示す。最後に、乱流への理論的アプローチとして、

ホップ発生点付近の振る舞いを近似する縮約方程式である複素ギンツブルグランダウ方程式を具体的に導出する。そして縮約方程式を利用して、拡散レプリケータ方程式の乱流解の統計則が特徴的なスケーリング則を持つことを示す。

2.1 モデルの説明と相図の作成

レプリケータ方程式についていくつか拡張したものが提案されたが、それらの一つとして種の移動を考慮に入れたモデルが提案された[97, 99, 100]. この試みは、拡散項をレプリケータ方程式につけることにより、詳細に議論された。拡散を伴うレプリケータ方程式は1989年にVickersにより提案された[99]。2種の場合、この方程式はフィッシャー方程式になる [101]。ここでフィッシャー方程式は侵略遺伝子の生成を記述する決定論モデルとして知られている [101, 102, 103, 104, 105]。ただし、3種以上の拡散レプリケータ方程式についてあまりに少ない情報しか得られてこなかった。(例外として2次元3種のレプリケータ方程式についての解析[106]と特殊な拡散の7種のレプリケータ方程式についての解析 [107]が挙げられる)。以下我々は今回、3種の拡散レプリケータ方程式を扱うことにする。

(24)

次の3種の拡散レプリケータ方程式を扱う、











∂X_i

∂t =X_i(

∑3 j=1

g_ijX_j −

∑3 j=1

∑3 k=1

g_jkX_jX_k) +D ∂²

∂r²X_i

∑3 i=1

X_i(r, t) = 1

(2.1)

(0≤X₁(r, t), X₂(r, t), X₃(r, t)≤1, [0≤r≤L])

ここでX_i =X_i(r, t)は種(又はプレイヤー)の頻度を表し(i= 1,2,3)、G={g_ij}は相互作用マトリクス、r ∈[0, L]は周期境界の1次元空間、Lはシステムサイズ、そして Dは拡散係数を表す。

簡単の為、相互作用行列G = {g_ij}について下のような非対称パラメータの場合に限定する。

G=







0 −1 1

1 0 −1

−1 1 0





+α







0 1 0 1 0 0 0 0 0





 (2.2)

ここでα は反対称（巡回）行列に加える対称摂動の強さを表す分岐パラメータである。

今回のパラメータについて説明すると、種X₁ < 種X₂ <種X₃ <種X₁・・・という力関係になるジャンケンゲームの力関係の均衡を、非対称パラメータαにより乱していくことになる。α= 0の時は、完全に対称なジャンケンゲームになっている。

他のパラメータDとLは固定する(D= 0.0009,L= 163.84)。つまりαは(2.1)式の唯一のパラメータとなる。なお今回のパラメータは、空間変換(r^′ =ar, a= 100/3 = 33.33· · ·) により、D = 1、L = 163.84/0.03 = 5461.33· · ·の場合と一致するということを述べておく。

最初に、(2.1)式の振る舞いを線形安定性解析を用いて予想してみる。(2.1)式には5 つの空間一様解が存在する。

(X₁⁰, X₂⁰, X₃⁰)₁ = (1,0,0), (2.3) (X₁⁰, X₂⁰, X₃⁰)2 = (0,1,0), (2.4) (X₁⁰, X₂⁰, X₃⁰)₃ = (0,0,1), (2.5)

(25)

図 1: 分岐ダイアグラム (α₁ = 0, α₂ ≃ −1.201) (X₁⁰, X₂⁰, X₃⁰)₄ = ( 1

3 +α, 1

3−α, 3−α²

(3 +α)(3−α)) (2.6) [√

3> α >−√ 3], (X₁⁰, X₂⁰, X₃⁰)₅ = (α−1

2α ,α+ 1

2α ,0) (2.7)

[α >1 又は −1> α].

(2.1)式への線形安定性解析の詳細はAppendix Aで扱っている．その結果を用いると、

大まかに次のことが予想される。我々の系はα≥√ 3、 √

3> α >0、 −1> αの3つのパラメータ領域でそれぞれ別の空間一様解が安定になる。しかし0≥ α ≥ −1の領域では空間一様安定解が一つも存在しない。したがってこの領域では我々の系は非定常な解をアトラクターに持つ。

大まかな予想だけでなく、詳細にわたる大域的な振る舞いや非線形効果も考慮に入れた振る舞いを調べる為、いくつかの差分法を用いてコンピュータシミュレーションを行う。シミュレーション手法の詳細はAppendix Bで述べる。そして数値計算の結果次のような解が現れた。α≥ α₁(= 0)、α < α2(∼=−1.201)の時、空間一様解が大域的に安定で一意に現れ、α1 ≥ α ≥ α₂では乱流解がアトラクターとして現れた。αについての相図は図1のようになる。乱流領域−1.19> α > α2では、初期値に応じていくつかの非乱流な局所解が大スケールの乱流アトラクタと共存して現れる。

α = −0.1、α = −1における乱流のスナップショットを図2 (a)(b)に載せる。そ

れらのスナップショットの3次元埋め込みプロット(ある時刻tでの、各rについての

(26)

図 2: X1(r, t)のスナップショット(tは固定); (a) はα=−0.1の場合で、(b)はα= −1 の場合である。 αを小さくすると振幅が大きくなる。これらのスナップショットの3次元空間への埋め込みプロット; (c)、(d)はそれぞれ、α= −0.1とα= −1の場合である。

(27)

(X1(r), X2(r), X3(r))の値をすべて3次元上に埋め込み、rについて隣同士の点を線でつないだ図)を図2 (c)(d)にのせる。ここで3次元埋め込みプロットは固定された時間 tにおけるX_i(r, t), (i= 1,2,3)のすべてのrの点から作られる。

それら数値シミュレーョンの結果は上記の線形安定性解析で得られた予測と概ね合っている。α=−1では、(2.3)式と(2.7)式の空間一様解が一致する。またα=−1付近では(2.7)式の不安定一様解が(2.3)式の安定一様解近くに残る。それゆえ、−1> α≥α₂ で乱流が現れ、(2.3)式の安定一様解に吸い込まれていかない。α =α₁(スーパークリティカルホップ分岐点)付近では、短い時間スケールでは安定振動が現れるように見えるが、長い時間スケールではそれが乱流解(α < α₁の時)か、または空間一様解（α > α1）に収束していく。これは、短い時間スケールではスーパークリティカルホップ分岐の影響が現れることを意味する。

2.2 乱流領域でのスペクトル解析

我々の主な関心は図1で現れる乱流領域(α1 > α > α2)の統計的性質の研究にある。

乱流の発生メカニズムについて、相関長と相関時間を用いた詳細な評価を行う。

2.2.1 空間スペクトル

システムサイズLを∆の間隔でNメッシュにわける、つまり∆ =L/N となる、そして空間スペクトルS_X_i(k)を以下のように定義する、

S_X_i(k) = 1 T

∫ _T

0

1

√N

∑N m=1

χ_i(∆m, t) exp (−imk)²dt, (2.8) ここでχ_i(r, t)はi種の密度X_i(r, t)の揺動部分であり(χ_i(r, t) =X_i(r, t)−< X_i(r, t)>_r、

< X_i(r, t) >_r= ¹_r∫L

0 X_i(r, t)dr)、X˜_i(k)はχ_i(r, t)の離散フーリエ変換、そして時間積分について実際には1000 ∼ 10000のサンプル時間から時間平均を行う。

乱流領域の両端近くの2つのケース（図1において）がそれぞれ図3 (a)、図3(b)である。分岐点付近(α₁−α ≤10⁻¹) では、パワースペクトルS_X₁(k)、SX2(k)、そして

(28)

SX3(k)はほとんど同じになる(図3(a))、

図3: 空間スペクトルS_X_i(k)の片対数プロット; (a)α=−0.1について、S_X₁(k)、SX2(k)、

そしてS_X₃(k)はほとんど同じ振る舞いを示す。(b)α=−1について、S_X₃(k)はS_X₁(k)、

S_X₂(k)とは異なる。これら二つの図より、波数が増える時、SXi(k)は「ほぼ」指数関数的に減衰する（単独の指数関数ではない）。ki1、ki2の詳細は(2.9)式で与えられる。

しかしパラメータαを減らす時（つまり、α₁−α≥10⁻¹の時）、S_X₃(k)は注目すべき違いを示す（図3(b)）。SXiの局所的な揺らぎにも関わらず、パワースペクトルS_X_i(k) が3つのタイプの特徴的な振る舞いを見せるのは、注目に値する。低波数領域(k < k_i1) では、スペクトルS_X_i(k)がストレッチドエクスポネンシャル形を取る。しかし高波数領域(k > k_i2)では、純粋な指数形を取る。そしてその間の領域(k_i1 < k < k_i2)では、

別の指数形を取る。つまり、

S_X_i(k) =











A_iexp(−d_ik^βⁱ) (k ≤k_i1) A^′_iexp(−d^′_i(k−k_i1)) (k_i1 < k≤k_i2) A^′′_i exp(−d^′′_i(k−k_i2)) (k_i2 < k).

(2.9)

となる。ここでk_i1 と k_i2 は図3で現れたクロスオーバー点であり、di, d^′_i, d^′′_i そして β_iはS_X_i(k)の数値計算の結果より得られたパラメータである(1.4 < β_i < 2)。同じスペクトルの振る舞いは、1次元cGL方程式で得られている、しかしそのスペクトルはストレッチドエクスポネンシャル的ではなく純粋に指数的である[108].

パラメータαを前から第2転移点(α≃α₂)に動かす時、(2.9)式で記述される3つの

(29)

図 4: (2.9)式の現象的係数のパラメータ依存性; (a) di, (b) d^′_i と d^′′_i, (c) βi, そして(d) ki1 と ki2.

特徴的領域は明確に観測されるが、現象的パラメータd_i, d^′_i, d^′′_i, β_i, k_i1, k_i2は比較的大きい揺らぎを見せる。さらに、図4で見られることだが、それぞれのパラメータが少なくとも第1転移点(α₁−α.10⁻¹)付近である種のスケーリング関係を満たすことは特筆に値する。

d_i ∝(α₁−α)⁻^0.786^±^0.001, (2.10)

d^′_i, d^′′_i ∝(α₁−α)−0.445±0.005, (2.11) k_i1, k_i2 ∝ |α₁−α|^0.44^±^0.01. (2.12) 図5は以下のように定義されるトータルパワーS_X_iを表す、

S_X_i ≡∑

k

S_X_i(k). (2.13)

ここで非-乱流領域(α > α₁, or α < α₂)ではトータルパワーは0である。なぜなら両方の領域でアトラクターは一様解であるからである。S_X_iは乱流の発生点付近(α≃α₁)で

(30)

図 5: (2.13)式で定義されるトータルパワーS_X_i

ほとんど線形である、しかしα ≃α₂付近で突然壊れる。SXiのαについての連続性と

(A.4)式の、分岐点で分散関係が純粋な虚数値になることより、α=α₁付近でスーパー

クリティカルホップ分岐が起きていると言える。α =α₂はサブクリティカル分岐で、

他のメタ安定アトラクターとカップルした逆分岐である崩壊型分岐が起きる。α =α₂ の分岐機構はまだ明らかになっていない。

2.2.2 相関長

前の結果を元に、3つの異なる方法を用いて乱流についての相関長(または特性長) の評価を行う。

1つ目は自己相関関数C_X_i(r)の半減長である。ここで、自己相関関数はスペクトル S_X_i(k)のフーリエ変換より、数値計算的に得られる（ウィナーヒンチンの定理[64]に依る）。図6(a)は規格化された相関関数C_X_i(r)/C_X_i(0)である。相関長r_iは

C_X_i(r_i) C_X_i(0) = 1

2 (2.14)

により求められる。もし図6(a)の相関減衰を単純な指数関数で近似すると、

C_X_i(r)∝exp(−r

ξ_i). (2.15)

Analysis of Heteroclinic Turbulence in Nonlinear Reaction Diﬀusion Equation

非線形反応拡散方程式における ヘテロクリニック乱流の解析