数値例におけるヘテロクリニック軌道の捕捉方法

6.2 数値例におけるヘテロクリニック軌道の捕捉方法

図9: _局所Lyapunov関数の定義域D_L11 と領域DA．DA上の点は必ずt → ∞^でx^∗_に収 束する．

6.2.2 x_u の捕捉方法

x_u ∈W^u(0)となる点x_u の捕捉方法について説明する．

(6.2)の定義域DLは第5.3節によれば

DL ={x| −0.1875≤x₁ ≤0.1875,−0.1875≤x₂ ≤0.1875}

である．ここで，Lyapunov関数値が0となる等位面 L₀ = {x ∈DL | L(x) = 0} ^を考え る．このとき，局所Lyapunov関数(6.2)からして

L₀= {x∈D_L |x₂ = 0}

だとわかる．さらに，Lyapunov関数値が正となる領域D_L⁺ = {x ∈ DL | L(x) > 0}^と Lyapunov関数値が負となる領域D⁻_L = {x∈DL |L(x)<0}^{を考える．このとき，}

D⁺_L ={x∈DL |x₂ <0}, D_L⁻ ={x∈DL |x₂ >0}

となる．平衡点以外でdL/dt < 0となることを考えれば，原点の不安定多様体W^u(0)上の 点x_u を捕捉するにあたって，領域D⁻_L と等位面L₀ に着目するのが妥当である．

次に，領域D_L⁻の境界∂D⁻_L :

∂D⁻_L =D_L1⁻ ∪D⁻_L2 ∪D_L3⁻ ∪L₀,

D_L1⁻ ={x∈D⁻_L |x₁ =−0.1875,0< x₂ ≤0.1875}, D_L2⁻ ={x∈D⁻_L | −0.1875< x1 <0.1875, x2 = 0.1875}, D_L3⁻ ={x∈D⁻_L |x1 = 0.1875,0< x2 ≤0.1875}

におけるフローについて考える．(6.1)からD_L1⁻ , D_L3⁻ 上の点ではx˙₁ > 0となり，D_L2⁻ 上 の点ではx˙₂ <0となる．L₀ においてはx˙₁ > 0,x˙₂ >0となる．局所Lyapunov関数の定

義域では閉軌道が存在しないことを鑑みれば，原点および点 (0.1875,0)^T を除くL0 上の点 から出発した軌道は正の有限時間でD_L⁻に流入し，D⁻_L3を通過して流出することがわかる．

図10: DL, D_L⁺, D⁻_L, D_L1⁻ , D⁻_L2, D_L3⁻ , L₀ の位置関係とD⁻_L1, D⁻_L2, D_L3⁻ , L₀ でのフローをま とめた図．黄色の矢印がD_L1⁻ , D⁻_L2, D⁻_L3 におけるフロー，赤色の矢印がL₀におけ るフローを表している．L₀\{0,(0.1875,0)^T}^上の点はD_L⁻ に流入し，D⁻_L3 上では D⁻_L から流失することがわかる．

以上の事実を用いて，次の手順でx_u の捕捉を行う．

1. 原点を挟んだL0 上の2点x_A,x_B :

x_A = (x^A₁,0)^T, x_B = (x^B₁ ,0)^T, x^A₁ <0< x^B₁

を選び，x_A,x_{B を初期点とする}(6.1)の解軌道ϕ(t,x_A), ϕ(t,x_B)を精度保証で算定 する．

2. ϕ(t,x_A), ϕ(t,x_B) がD_L3⁻ を通過する時刻 ti, Tj > 0 を探し出す．これは， ODE-IVPの精度保証法を用いて得られるϕ(t,x_A), ϕ(t,x_B)の時刻ti, Tj における解区間 [ϕ(t_i,x_A)],[ϕ(T_j,x_B)] :

[ϕ(ti,x_A)] :={x∈R² | x_i ≤x₁ ≤xi, y_i ≤x₂ ≤y_i}, [ϕ(Tj,x_B)] :={x∈R² |x_j ≤x1≤ xj, y_j ≤x2 ≤y_j},

x_i, xi, y_i, y_i, x_j, xj, y_j, y_j ∈R について，

x_i, x_j >0.1875 (6.5)

を示せればよい．

3. ϕ(t,x_A), ϕ(t,x_B)がD_L3⁻ を通過する前の時刻ti−1, Tj−1 > 0を探し出す．これは，

ODE-IVPの精度保証法を用いて得られるϕ(t,x_A), ϕ(t,x_B)の時刻t_i−1, T_j−1にお ける解区間[ϕ(t_i−1,x_A)],[ϕ(T_j−1,x_B)] :

[ϕ(t_i−1,x_A)] :={x∈R² |x_i−1 ≤x₁ ≤x_i−1, y_i−1 ≤x₂ ≤y_i−1},

[ϕ(Tj−1,x_B)] :={x∈R² |x_j−1 ≤x1 ≤xj−1, y_j−1 ≤x2 ≤y_j−1}, x_i−1, x_i−1, y

i−1, y_i−1, x_j−1, x_j−1, y

j−1, y_j−1 ∈R について，

x_i−1, x_j−1 <0.1875 (6.6) を示せればよい．

図11: 左図は解区間[ϕ(Ti,x_A)],[ϕ(t_i−1,x_A)]とD_L⁻, D⁻_L3 との位置関係を,右図は 解区間 [ϕ(Tj,x_B)],[ϕ(T_j−1,x_B)]とD_L⁻, D⁻_L3との位置関係を描いている．

4. ϕ(t,x_A)について，適当な時間分点0 = t₀ < t₁ < · · · < t_i−1 < ti を設定し，時刻 [0, t1],[t1, t2],· · ·,[ti−1, ti]における解区間

[ϕ([t₀, t₁],x_A)],[ϕ([t₁, t₂],x_A)],· · · ,[ϕ([t_i−1, ti],x_A)] (6.7) を精度保証で求める．また，ϕ(t,x_B) についても適当な時間分点 0 = T₀ < T₁ <

· · ·< T_j−1 < Tj を設定し，時刻[0, T₁],[T₁, T₂],· · · ,[T_j−1, Tj]における解区間 [ϕ([T₀, T₁],x_B)],[ϕ([T₁, T₂],x_B)],· · ·,[ϕ([T_j−1, T_j],x_B)] (6.8) を精度保証で求める．区間時間における解区間は，第2.4.4節で記した Taylor展開 による精度保証法(2.5)におけるステップ幅hを区間時間に置き換えることで得られ る．なお，kvライブラリ([11])にはこの機能が実装されている．

このとき，

[ϕ(t_k−1,x_A)],[ϕ(tk,x_A)]⊂[ϕ([t_k−1, tk],x_A)],

[ϕ(T_l−1,x_B)],[ϕ(Tl,x_B)]⊂[ϕ([T_l−1, Tl],x_B)], (6.9) 1≤k ≤i, 0≤l ≤j

であり，また

[ϕ([t_k−1, tk],x_A)]∩[ϕ([tk, t_k+1],x_A)]6=∅,

[ϕ([T_l−1, Tl],x_B)]∩[ϕ([Tl, T_l+1],x_B)]6=∅, (6.10) 1≤k ≤i−1, 0≤l ≤j −1

となることに注意されたい．

5. 得られた解区間(6.7),(6.8)をもとに，解区間の和集合A, B :

A= [ϕ([0, t₁],x_A)]∪[ϕ([t₁, t₂],x_A)]∪ · · · ∪[ϕ([t_i−1, ti],x_A)], B= [ϕ([0, T₁],x_B)]∪[ϕ([T₁, T₂],x_B)]∪ · · · ∪[ϕ([T_j−1, Tj],x_B)]

を考える．このとき，Aはϕ(t,x_A)の時刻0 ≤t ≤ti における値の存在範囲を，B はϕ(t,x_B)の時刻0 ≤t≤ Tj における値の存在範囲を意味する．また，(6.10)より A, B は連結集合である．ここで，次の2つの条件を満たしているか確認する．

A\{[ϕ([ti−1, ti],x_A)]}, B\{[ϕ([Ti−1, Ti],x_A)]} ⊂(D⁻_L ∪L0) (6.11)

 A∪B=∅ (6.12)

(6.11) については，ϕ(t,x_A), ϕ(t,x_B) が時刻 0 ≤ t ≤ t_i−1,0 ≤ t ≤ T_j−1 で局所

Lyapunov関数(6.2)の定義域に含まれているかを確認するためである．解析的には

領域D⁻_L からはD⁻_L3以外の境界で流出しないことがわかっているが，区間演算にお ける誤差拡大によって条件(6.11)を満たさない場合が存在し得る．

A, B はともに連結集合であることや条件 (6.11),(6.12)から，集合A, B の位置は以下の図 のように描ける．

図12: _条件(6.11),(6.12)を満たすときのDL, D⁻_L, D_L3⁻ , L₀ と点x_A,x_{B と集合}A, Bの位 置関係．黄色の領域が集合A,緑色の領域が集合Bである．

解区間[ϕ([t_i−1, ti],x_A)],[ϕ([T_j−1, Tj],x_B)]を

[ϕ([t_i−1, ti],x_A)] :={x∈R² | x_ti ≤x₁ ≤xti, y

ti ≤x₂ ≤y_ti} [ϕ([T_j−1, Tj],x_B)] :={x∈R² | x_{T j} ≤x₁ ≤xT j, y_{T j} ≤x₂ ≤y_{T j}},

x_{T i}, xT i, y

T i, y_{T i}, x_{T j}, xT j, y

T j, y_{T j} ∈R

と表す．このとき，解区間 [ϕ([t_i−1, ti],x_A)],[ϕ([T_j−1, Tj],x_B)] と D_L3⁻ との共通部分 XLA,XLBは

XLA= [ϕ([t_i−1, ti],x_A)]∪D_L3⁻

={x∈R² | x₁ = 0.1875, y

ti ≤x₂ ≤y_ti}, XLB = [ϕ([Tj−1, Tj],x_B)]∪D_L3⁻

={x∈R² | x1 = 0.1875, y_{T j} ≤x2 ≤y_{T j}} となる．さらに，X_LA,X_LB の和集合を次の区間ベクトル[X_u] :

[Xu] ={x∈R² | x₁ = 0.1875, y_{T j} ≤x₂ ≤y_ti}

= ([0.1875,0.1875],[y_{T j}, y_ti])^T

を考える．このとき，図12,局所Lyapunov関数の定義域(6.3)では閉軌道が存在しないこ と(付録B.3を参照)，自励系において解軌道は交叉しないこと(付録C.1を参照)から，

∃x_u ∈[X_u], lim

t→−∞x_u =0 が成立する．また，[Xu]がx_u の存在範囲となる．

6.2.3 [Xu]からDAへの解軌道の検証

前節で得られた[Xu]を初期区間とする解軌道ϕ(t,[Xu])について，

[ϕ(TU S,[Xu])]⊂DA (6.13)

を満たす時刻 TU S > 0を探し出す．なお，[ϕ(TU S,[Xu])]はODE-IVPの精度保証法を用 いて得られた，時刻TU S におけるϕ(TU S,[Xu])の解区間を表す．(6.13)の判定方法につい ては，解区間[ϕ(TU S,[Xu])]について

supL₁₁([ϕ(TU S,[Xu])])≤α

となればよい．なお，supL₁₁([ϕ(T_{U S},[X_u])])とはL₁₁([ϕ(T_{U S},[X_u])]の上端を表す．

このとき，[Xu]⊂W^s(x^∗)となる．これで，

∃x_u ∈[xu], lim

t→∞ϕ(t,x_u) =x^∗∧ lim

t→−∞ϕ(t,x_u) =0 が成立し，ϕ(t,[Xu])はヘテロクリニック軌道の存在範囲となる．