カオス・フラクタル 講義ノート #8
担当 : 井上 純一 (情報科学研究科棟 8-13)
URL : http://chaosweb.complex.eng.hokudai.ac.jp/˜j inoue/index.html

平成 23 年 6 月 21 日

目次
11 非線形力学系とカオス
11.1 レスラー方程式とその振る舞い . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70
71

11.2 ポアンカレ断面 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

課題 6 の解答例

dv
dt
dxdt

= (1 − x2 )v − x(167)= v(168)の振る舞いを数値的に調べてみる. まず, パラメータ  がゼロの場合, 方程式は単振動 (調和振動子)を表す. 横軸を x, 縦軸を v とした場合の (位相空間内の) 軌道は (167) 式の両辺に v = dx/dt をかけるとエネルギーの時間変化についての関係式がでるので, 実行するとdxdt(169)( 2)v
x2+=022(170)v

dv
dt= −xすなわちddtとなり, E を初期条件から決まる一定値とすると

v 2 + x2=√( 2E)2(171)√2E = v02 + x20 の円となる (x0 , v0 は t = 0 のときの位置と速度). 実際に  = 0の場合に (167)(168) 式を 4 次のルンゲ・クッタ法で数値的に解き, 軌道 (x, v) をプロットすると図となり, 半径が√39 のようになる. これを計算するためのソースコードを次に載せる.
#include <stdio.h>#include <math.h>#include <stdlib.h>#define h 0.01
68

カオス・フラクタル 2011 #8

担当: 大学院情報科学研究科 井上 純一

4ε = 0: Harmonic

32

1

v
0-1-2-3-4-4-3-2

√

-1

0x

1

2

34

√

√
√
v02 + x20 = 12 + 32 = 10 = 3.16 とした円になる. gnuplot によってプロットした図は縦軸と横軸のスケールが違っていることに注意.

図 39:  = 0 の場合の軌道. 半径を 2E =#define epsilon 0.0
double func1(x,y)
double x;
double y;{return (y);}double func2(x,y)
double x;
double y;{return (epsilon*(1.0-x*x)*y-x);}main(){FILE *fpr;
int i,imax=50000;
double x,k1x,k2x,k3x,k4x,kx,y,k1y,k2y,k3y,k4y,ky,z,k1z,k2z,k3z,k4z,kz;
if((fpr = fopen("test0.dat", "wt")) !=NULL){
for(i = 0,x=1.01,y=3.01; i <= imax; i++){k1x = h*func1(x,y);k1y = h*func2(x,y);k2x = h*func1(x+0.5*k1x,y+0.5*k1y);
k2y = h*func2(x+0.5*k1x,y+0.5*k1y);
k3x = h*func1(x+0.5*k2x,y+0.5*k2y);
k3y = h*func2(x+0.5*k2x,y+0.5*k2y);
k4x = h*func1(x+k3x,y+k3y);
k4y = h*func2(x+k3x,y+k3y);
kx = (k1x+2.0*k2x+2.0*k3x+k4x)/6.0;

ここは 69 ページ目

カオス・フラクタル 2011 #8

担当: 大学院情報科学研究科 井上 純一

ky = (k1y+2.0*k2y+2.0*k3y+k4y)/6.0;
x = x + kx;
y = y + ky;
fprintf(fpr,"%lf %lf %lf\n",i*h,x,y);}}
fclose(fpr);}
いくつかの有限の  に対する位相空間内の軌道を求めてみると図のようになる. 図 40 (上, および43

ε = 0.01

ε = 0.5

32211v

v
00-1-1-2-2-3-4-4-3-2-10x

1243-3-2.5

4-2-1.5-1-0.5

0x 0.511.5

3

ε = 1.5

22.5x, ε = 0.5v, ε = 0.532211v
00-1-1-2-2-3-4-2.5-2-1.5-1-0.5

0x 0.5

11.5

22.5-3050

100150t200

図 40: 有限の  の場合のアトラクタの様子と  = 0.5 の場合の時間変化 (右).
左下) は  = 0.01, および  = 0.5, 1.5 の場合の軌道をそれぞれプロットした. また, 右下は  = 0.5の場合の位置 x, 速度 v の時間変化である. これらの図より, 時間の経過とともに, 軌道はある閉曲線に収束することがわかる. このような軌道の収束する点, もしくは曲線をアトラクタと呼ぶ.

11

非線形力学系とカオス

今回の講義では, 前回の講義で学んだレスラー方程式の振る舞いを方程式の形から定性的に説明するとともに, アトラクタを可視化するためのポアンカレ断面について詳しくみていく.

ここは 70 ページ目

カオス・フラクタル 2011 #8

11.1

担当: 大学院情報科学研究科 井上 純一

レスラー方程式とその振る舞い

前回みたように, レスラー方程式:

dx
dt
dy
dtdz
dt

= −y − z(172)= x + ay(173)

= b + z(x − c)

(174)

をルンゲ・クッタ法を用いて数値的に解き, その軌道を 3 次元に描画すると図 41 のようになる. こa=0.398, b=2, c=465z43210-4-3-2-10x123456 -6-5-4-3-2-10123y

図 41: レスラー方程式の数値解から得られる軌道の様子. パラメータは a = 0.398, b = 2, c = 4 に選んである.の図を描画するのに, 方程式の初期値は原点近傍 (x, y, z) = (0.01, 0.01, 0.01) に選んである. この図をみるからに複雑なアトラクタであることがわかるが, この振る舞いを方程式の形からある程度説明することができる.
まず, この図 41 から, 軌道がスタートしてからしばらくは x-y 平面内に停滞していることがわかる. そこで, 方程式 (172)(173)(174) において, x-y 平面内近傍での初期の振舞いがどうであるかを調べるために, z ' 0 とおいてみると, (172)(173) 式から

dxdtdydt= −y(175)= x + ay(176)

ここは 71 ページ目

カオス・フラクタル 2011 #8

担当: 大学院情報科学研究科 井上 純一が得られる. (175) 式の両辺をを時間で微分し, それに (176) 式を代入すると, 次の x に関する 2 階の線形微分方程式が得られる.

d2 xdx−a+x = 0dt2dt(177)

そこで, これに x = ept を代入すると, 特性方程式は p2 − ap + 1 = 0 となるので, この 2 次方程式の解はp =であり, 0 < a < 2 のときに√a±√a2 − 42(178)

4 − a2 /2 ≡ ω とおけば, この微分方程式 (177) の解は{}
x = eat/2 αeiωt + βe−iωt = eat/2 {A cos ωt + B sin ωt}(179)となる. 従って, y は (176) 式から

ydx= −= −eat/2dt{()()}
aAaB+ Bω cos ωt +− Aω sin ωt
22(180)である (A, B は初期条件により決まる積分定数). 従って, x-y 平面内での軌道はその「半径」が時間とともに eat/2 で大きくなって行くような円運動 (らせん運動) を行う.
ここで, z の寄与が何も無ければ, らせん運動の半径は時間の経過とともに指数関数的な速さで無限大に達するはずであるが, 図 41 からわかるように, そうはなっていない. 従って, 時間の経過とともに z からの寄与が効いてきているはずなので, それを調べるため, (174) 式を見てみる. この式より, x の値が c 以下であれば, dz/dt = 0 となる点: z = −b/(x − c) > 0 へ収束する. しかし, 既にみたように, x の絶対値は時間とともに大きくなるのであるから, らせん軌道上で x > c となる
点で dz/dt = b + z(x − c) > 0 となり (b > 0 と選んでいるので (ちなみに図 41 では b = 2)), z はこの点で急激な増加に転じる (この増加の「急激さ」は b が大きいほど大きい). 実際に図 42 を見ると, 確かに x の大きさがある値を超えた瞬間に z 方向への軌道の急激な増加が見られる.
しかし, 同時にこの図 42 より, この z の方向の増加は長続きせず, 軌道は x-y 平面に引き戻されている. これは方程式からも見て取れる. つまり, z が大きくなると (172) 式から dx/dt の符号が負に傾き始め, 従って, x が減少し始める. すると, z が増加するための条件であった x > c が崩れ始め, やがて x < c となり, dz/dt = b + z(x − c) < 0 へと転じることから z が減少しはじめ, やがて, x-y 平面に引きつけられることになる. 上記の一連の挙動を繰り返すことで図 41 が得られることになる.

11.2

ポアンカレ断面

前節では数値実験によって得られた「奇妙な」アトラクタ (実際に奇妙なアトラクタ, あるいはストレンジ・アトラクタと呼ばれる) が, その方程式の振る舞いからある程度までは定性的に説明できることを見た. 図 41 に示した奇妙なアトラクタは 3 次元空間内の軌道であるので, 角度を変えてみることを繰り返すことにより, その構造を視覚的に観察することが重要である. 例えば, gnuplot
ではgnuplot> set view 60,15

ここは 72 ページ目

カオス・フラクタル 2011 #8

担当: 大学院情報科学研究科 井上 純一

z43.5
32.5
21.5
10.50-3-2-1012x345 -5-4-3-2-1012y

図 42: 図 41 の軌道のうち, ルンゲ・クッタ法でのステップを途中で打ち切って z が急激に増加する部分だけを抜き出したもの. x > c となる時点で z 方向への急激な増加が始まる.

で見る角度を指定することができる (60,15 などの数値を適時変える. これは後の演習の時間に詳しく説明する).
しかし, ここでは, 軌道をある平面で切った切り口を考えてみる. 3 次元軌道を平面で切ったの断面なので, 得られるものは 2 次元像である. そこで, x-y 面内における x 軸からのなす角度 θ を用いて y cos θ − x sin θ = 0 とし, この平面と z 軸で囲まれた「半平面」でこの軌道を切ってみる. その結果を図 43 に載せよう. この図 43 (図 41 と同じ図) に示した軌道を y cos θ − x sin θ = 0 と z 軸で囲まれる半平面で切った断面が右図である. この図より, 角度 θ を θ = 0 から θ = 300o へ増やすにつれて, 軌道の断面を表す「線分」が「折りたたまれ」, 次いで「引き伸ばされる」様子が見て取れる.
このように, 奇妙なアトラクタをある半平面で切り取った断面をポアンカレ断面と呼ぶ. 後の演習の時間で実際に求めてもらうが, このポアンカレ断面を通過する点列を (xt , xt+1 ) としてプロットすることで, この力学系を表す実質的な「写像」を取り出すことができる1 .
さて, 上に述べた奇妙なアトラクタを示す力学系 (微分方程式) はレスラー方程式だけではない.
例えば, 次のようなローレンツ方程式:

dx
dt
dy
dtdz
dt

=−ax + ay(181)=cx − y − xz(182)=−bz + xy(183)

もパラメータ a, b, c の選び方によっては奇妙なアトラクタを持つ. その結果の一例を図 44 に載せる. 講義ではこのローレンツアトラクタの動画デモを見せる予定である.
1 「リターンマップ」と呼ばれる.

ここは 73 ページ目

カオス・フラクタル 2011 #8

担当: 大学院情報科学研究科 井上 純一

a=0.398, b=2, c=4654z

32
10-4-3-2-10

1x23456 -6-5-4-3-2

4.5-1

0123

y

θ=0
θ = 60θ = 120
θ = 180
θ = 240
θ = 300

43.53z

2.5
2
1.5
10.5
00

123r

4567

図 43: 左図の軌道を y cos θ − x sin θ = 0 と z 軸で囲まれる半平面で切った断面が右図.

レポート課題 8ローレンツ方程式を a = 10, b = 8/3, c = 30 に対して数値的に解き, 平面 z = 20, 30, 40 で切り取ったポアンカレ断面を求めるプログラムを作成し, 描画ソフトが使えるものはそれをプロットせよ (その手のソフトが使えない者は後の演習の時間に説明するので, 現段階ではソースコードのみで良い).

ここは 74 ページ目

カオス・フラクタル 2011 #8担当: 大学院情報科学研究科 井上 純一

a=10, b=8/3, c=28z50
4540353025201510
50-20-15-10-5x0

51015-10-15-20
20-25-50

51015202530y

図 44: ローレンツ方程式 (181)(182)(183) から得られる奇妙なアトラクタ. a = 10, b = 8/3, c = 28 に選んである.

ここは 75 ページ目