とともにzからの寄与が効いてきているはずなので,それを調べるため, (184)式を見てみる. この 式より,xの値がc以下であれば,dz/dt= 0となる点: z=−b/(x−c)>0へ収束する. しかし,既 にみたように, xの絶対値は時間とともに大きくなるのであるから, らせん軌道上でx > cとなる 点でdz/dt=b+z(x−c)>0となり(b >0と選んでいるので(ちなみに図44ではb = 2)),zは この点で急激な増加に転じる(この増加の「急激さ」はbが大きいほど大きい). 実際に図45を見
-3 -2 -1 0 1 2 3 4 5-5
-4 -3 -2 -1 0 1 2
0 0.5 1 1.5 2 2.5 3 3.5 4 z
x
y z
図 45: 図44の軌道のうち,ルンゲ・クッタ法でのステップを途中で打ち切ってzが急激に増加する部分だけを抜き出し たもの.x > cとなる時点でz方向への急激な増加が始まる.
ると,確かにxの大きさがある値を超えた瞬間にz方向への軌道の急激な増加が見られる.
しかし,同時にこの図45より,このzの方向の増加は長続きせず,軌道はx-y平面に引き戻され ている. これは方程式からも見て取れる. つまり,zが大きくなると(182)式からdx/dtの符号が 負に傾き始め,従って, xが減少し始める. すると, zが増加するための条件であったx > cが崩れ 始め,やがてx < cとなり, dz/dt=b+z(x−c)<0へと転じることからzが減少しはじめ,やが て,x-y平面に引きつけられることになる. 上記の一連の挙動を繰り返すことで図44が得られるこ とになる.
で見る角度を指定することができる(60,15などの数値を適時変える. これは後の演習の時間に詳 しく説明する).
しかし,ここでは,軌道をある平面で切った切り口を考えてみる. 3次元軌道を平面で切ったの断 面なので,得られるものは2次元像である. そこで,x-y面内におけるx軸からのなす角度θを用い てycosθ−xsinθ= 0とし,この平面とz軸で囲まれた「半平面」でこの軌道を切ってみる. その 結果を図46に載せよう. この図46(図44と同じ図)に示した軌道をycosθ−xsinθ= 0とz軸で
-4 -3 -2 -1 0 1 2 3 4 5 6-6-5-4-3-2-1 0 1 2 3 0
1 2 3 4 5 6 z
a=0.398, b=2, c=4
x
y z
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
0 1 2 3 4 5 6 7
z
r
θ = 0 θ = 60 θ = 120 θ = 180 θ = 240 θ = 300
図46: 左図の軌道をycosθ−xsinθ= 0とz軸で囲まれる半平面で切った断面が右図.
囲まれる半平面で切った断面が右図である. この図より,角度θをθ= 0からθ= 300oへ増やすに つれて,軌道の断面を表す「線分」が「折りたたまれ」, 次いで「引き伸ばされる」様子が見て取 れる.
このように,奇妙なアトラクタをある半平面で切り取った断面をポアンカレ断面と呼ぶ. 後の演 習の時間で実際に求めてもらうが,このポアンカレ断面を通過する点列を(xt, xt+1)としてプロッ トすることで,この力学系を表す実質的な「写像」を取り出すことができる35.
さて,上に述べた奇妙なアトラクタを示す力学系(微分方程式)はレスラー方程式だけではない.
例えば,次のようなローレンツ方程式:
dx
dt = −ax+ay (191)
dy
dt = cx−y−xz (192)
35「リターンマップ」と呼ばれる.
dz
dt = −bz+xy (193)
もパラメータa, b, cの選び方によっては奇妙なアトラクタを持つ. その結果の一例を図47に載せ る. 講義ではこのローレンツアトラクタの動画デモを見せる予定である.
-20 -15 -10 -5 0 5 10 15 20-25-20-15-10-5 0 5 10 15 20 25 30 0 5
10 15 20 25 30 35 40 45 50
z
a=10, b=8/3, c=28
x
y z
図47: ローレンツ方程式(191)(192)(193)から得られる奇妙なアトラクタ. a= 10, b= 8/3, c= 28に選んである.
レポート課題8
ローレンツ方程式をa= 10, b= 8/3, c= 30に対して数値的に解き,平面z= 20,30,40で切り取っ たポアンカレ断面を求めるプログラムを作成し,描画ソフトが使えるものはそれをプロットせよ(そ の手のソフトが使えない者は後の演習の時間に説明するので,現段階ではソースコードのみで良い).
ここは ページ目
課題8の解答例
先週の復習として, レポート課題を確認しておこう. a= 10, b= 8/3, c= 30 と選んだ場合のロー レンツ方程式のz= 20,30,40におけるポアンカレ断面を図48に載せる. また,このグラフを作成
-20 -15 -10 -5 0 5 10 15 20
-20 -15 -10 -5 0 5 10 15
y
x
z = 20 z = 30 z = 40
図48: a= 10, b= 8/3, c= 30と選んだ場合のローレンツ方程式のz= 20,30,40におけるポアンカレ断面.
するには
x = x + kx;
y = y + ky;
z = z + kz;
arg = fabs(z-20);
if((i>=iSkip) && (arg <5.0e-6)){
fprintf(fpr,"%lf %lf\n",x,y);
}
のように,z= 20をある精度で満たす(x, y)のみをファイルに書き込むようにすればよい.
12 アトラクタの埋め込み次元と相関次元
前回までの講義ではレスラー方程式やローレンツ方程式などの非線形微分方程式の数値解が奇妙 なアトラクタを持つことをみたが,レスラー方程式はx-y平面内に強く引きつけられることで3次 元というよりももしろ2次元的な構造を持っている. そこで,ここではこの奇妙なアトラクタの実 質的な次元を評価するためのいくつかの方法をみていく.
12.1 埋め込み次元とアトラクタの再構成
レスラー方程式:
dx
dt = −y−z (194)
dy
dt = x+ay (195)
dz
dt = b+z(x−c) (196)
を数値的に解いた後,例えばx(t)のみに着目し,適当な時間遅れτを導入し,次のようなm次元ベ クトルを作る.
(x(t), x(t+τ),· · ·, x(t+ (m−1)τ)) (197) そこで,このm次元ベクトルを逐次プロットすることで,前回みた奇妙なアトラクタが再構成できるか を調べてみよう. 図49に次元をm= 2と選んだ場合の(x(t), x(t+τ)),(y(t), y(t+τ)),(z(t), z(t+τ)) をぞれぞれプロットした. この図より,xのみ, y,zのみの軌道を適切な時間遅れτ, 埋め込み次元
-4 -3 -2 -1 0 1 2 3 4 5 6
-4 -3 -2 -1 0 1 2 3 4 5 6
x(t+τ)
x(t)
τ = 2
-6 -5 -4 -3 -2 -1 0 1 2 3
-6 -5 -4 -3 -2 -1 0 1 2 3
y(t+τ)
y(t)
τ = 2
0 1 2 3 4 5 6
0 1 2 3 4 5 6
z(t+τ)
z(t)
τ = 1
図49: 再構成.埋め込み次元をm= 2に選んでいる.
mで切り出すと,x, y, zを3次元にプロットした場合とほぼ同じアトラクタを再構成することがで きる. この例ではm= 2,τ= 1,2に選んでいるが,与えられたデータが方程式からの時系列ではな く性質が不明なものである場合,これら2つのパラメータを適切に決定しなければならない. 時間 遅れτはx(t)とx(t+τ)の自己相関関数,あるいは,相互情報量を評価することで見積もることが できるが,ここではそれには触れない. 次元mに関しては次に述べる「相関次元」がその指針を与
ここは ページ目
える. そこで今回の講義では,相関次元の計算を見ていく. その準備として,まずは「次元」につい て確認しておこう.