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

平成 23 年 6 月 28 日

目次
12 アトラクタの埋め込み次元と相関次元12.1 埋め込み次元とアトラクタの再構成 . . . . . . . . . . . . . . . . . . . . . . . . . .

77
77

12.2 次元再考 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12.3 相関次元 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

787912.4 相関次元の一般化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83

課題 8 の解答例
先週の復習として, レポート課題を確認しておこう. a = 10, b = 8/3, c = 30 と選んだ場合のローレンツ方程式の z = 20, 30, 40 におけるポアンカレ断面を図 45 に載せる. また, このグラフを作成
20z = 20z = 30z = 40

15105y0-5-10-15-20-20-15-10-5x051015

図 45: a = 10, b = 8/3, c = 30 と選んだ場合のローレンツ方程式の z = 20, 30, 40 におけるポアンカレ断面.
するには76カオス・フラクタル 2011 #9

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

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

埋め込み次元とアトラクタの再構成

レスラー方程式:

dxdt
dy
dtdz
dt

= −y − z(184)= x + ay(185)

= b + z(x − c)(186)を数値的に解いた後, 例えば x(t) のみに着目し, 適当な時間遅れ τ を導入し, 次のような m 次元ベクトルを作る.

(x(t), x(t + τ ), · · · , x(t + (m − 1)τ ))(187)

そこで, この m 次元ベクトルを逐次プロットすることで, 前回みた奇妙なアトラクタが再構成できるかを調べてみよう. 図 46 に次元を m = 2 と選んだ場合の (x(t), x(t+τ )), (y(t), y(t+τ )), (z(t), z(t+τ ))をぞれぞれプロットした. この図より, x のみ, y, z のみの軌道を適切な時間遅れ τ , 埋め込み次元m で切り出すと, x, y, z を 3 次元にプロットした場合とほぼ同じアトラクタを再構成することができる. この例では m = 2, τ = 1, 2 に選んでいるが, 与えられたデータが方程式からの時系列ではなく性質が不明なものである場合, これら 2 つのパラメータを適切に決定しなければならない. 時間遅れ τ は x(t) と x(t + τ ) の自己相関関数, あるいは, 相互情報量を評価することで見積もることができるが, ここではそれには触れない. 次元 m に関しては次に述べる「相関次元」がその指針を与える. そこで今回の講義では, 相関次元の計算を見ていく. その準備として, まずは「次元」について確認しておこう.

ここは 77 ページ目

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

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

63

τ=25

τ=224

1

30-1y(t+τ)

x(t+τ)

21-20-3-1-4-2-5-3-4-4-3-2-101x(t) 2345-66-6-5-46-3-2y(t) -10123

τ=15z(t+τ)

432100123z(t)456

図 46: 再構成. 埋め込み次元を m = 2 に選んでいる.

12.2次元再考

空間の次元に関し, 我々は直線が 1 次元, 平面が 2 次元等を知っている. しかし, ここでは任意の対象の次元を少し変わった方法で定義する. この方法は後に「フラクタル」を学ぶ際にも重要である.
まず, 半径が  の球 (円) で問題とする対象物を覆うことを考えよう. このとき, 対象を隙間無く覆うために必要な球の個数を n() とおく. この n() を用いて, この対象物の次元を次式で定義する.

log n()→0 log = − lim

D0(188)この次元 — 容量次元と呼ばれる— が我々が経験的に知っている次元に反しないかを確認するた
め, 長さ l の線分の次元 D0 を計算してみよう.
まず, 長さ l の線分を隙間無く半径  の球で覆うためには l/2 個の球が必要であるから, n() = l/2である. 従って

D2 (長さ l の線分)log(l/2)= − lim→0→0log = − lim{log(l/2)
log }+log =1log (189)となり, 線分の次元は 1 であるという事実に合致する.
また, 1 辺の長さが l の正方形に対しては同様にしてこれを覆うための球の個数は n() = l2 /π2であるから

D2 (一辺の長さが l の正方形) =log(l2 /π2 )= − lim
→0→0log − lim{log(l2 /π)log }+2log = (2190)log 

ここは 78 ページ目

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

担当: 大学院情報科学研究科 井上 純一となり, やはり我々の経験事実に合う. 従って, この方法を用いることで任意の図形の次元を計算することができる. 既にみた奇妙なアトラクタの次元もこの考え方を拡張することで算出することができる.
ここで述べた 2 つの例はほとんど自明な場合であったのだが, 海岸線など複雑に入り組んだ図形の次元を上記の方法で見積もると 1, 2, 3, · · · のような整数次元ではなく, 非整数の次元を持つことがわかっている. これについては再来週からの [フラクタル編] で詳しくみていくことになる.

12.3

相関次元

前節で見たように, ある力学系から生成されたデータ列 {x(ti )}: i = 1, 2, · · · , N を時間遅れ τ でサンプリングし, 次の m 次元のベクトル:ri≡(x(ti ), x(ti + τ ), · · · , x(ti + (m − 1)τ )), i = 1, · · · , N(191)

を考える. このとき, r i を中心とする半径  の球を考え, この球内にアトラクタ上の点 r j : j =1, · · · , N が入る確率 pi を考えると
1 ∑Θ( − |r i − r j |)N j=1
Npi=(192)となる. ここで, Θ(x) は階段関数 (ステップ関数) と呼ばれる関数であり, x ≥ 0 で 1, x < 0 でゼロを返す. また, |r i − r j | を陽な形で書き出しておくと|r i − r j | = |x(ti ) − x(tj )| + |x(ti + τ ) − x(tj + τ )| + · · · + |x(ti + (m − 1)τ ) − x(tj + (m − 1)τ )|(193)であることに注意されたい. このとき相関次元を∑N
D2=lim

→02i=1 pilog (194)で定義する. ここでN∑
i=1p2i=N N
1 ∑∑Θ( − |r i − r j |) ≡ C()N 2 i=1 j=1(195)に注意すると相関次元 (194) は C() を用いて次のように書き直せる.

D2=log C()→0log lim

(196)

よって,  に関する極限をひとまず無視すれば, D2 = log C()/log , すなわちC()= D2(197)が成り立つ. 従って (195) 式に従って  の変化に対する C() を算出し, 得られる結果を両対数プロットした際のグラフの傾きが相関次元 D2 を与える.
レスラー方程式において x(t) から時間遅れ τ = 0.01 × 200 = 2 のもとで m = 1, 2, 3, 4, 5 次元のベクトル (x(ti ), x(ti + τ ), · · · , x(ti + (m − 1)τ )) : i = 1, · · · , N を作成し, それから上記レシピに
ここは 79 ページ目

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

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

0-1log C(ε)

-2

m=1
m=2
m=3
m=4
m=5-3-4-51

1.2

1.4log ε1.61.82

図 47: レスラー方程式において, x(t) から時間遅れ τ = 2 のもとで作成した m = 1, 2, 3, 4, 5 次元のベクトル (x(ti ), x(ti +τ ), · · · , x(ti + (m − 1)τ )) : i = 1, · · · , N から算出した相関次元. N = 3000 として計算してある. m = 2 とすると傾きD2 ∼ 2 の直線が得られ, この傾きは以降 m を増やしても変わらない.

従って相関次元を求めた結果を図 47 載せよう. この図より, m を増加するにつれ, グラフの傾きは実質的なアトラクタの次元である D2 = 2 に近づくことがわかる. この方法はレスラー方程式からの軌道にみならず, 何らかの方法で実測されたデータ, もしくは背後の力学系や方程式が未知である経済/金融時系列データなどにも適用可能な方法である.
参考までに m = 5 の場合のプログラムのコード例を載せておこう.

/********************************************************************//*
相関次元の計算：レスラー方程式*//********************************************************************/#include <stdio.h>#include <math.h>#define h 0.01
#define a 0.398#define b 2#define c 4#define imax 50000#define N 3000
double g1[N];double g2[N];double g3[N];double g4[N];double g5[N];double xx[imax];
double func1(x,y,z)

ここは 80 ページ目カオス・フラクタル 2011 #9

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

double x;
double y;
double z;{return (-(y+z));}double func2(x,y,z)
double x;
double y;
double z;{return (x+a*y);}double func3(x,y,z)
double x;
double y;
double z;{return (b+z*(x-c));}double STEP(x)double x;{if(x>=0.0){
return (1.0);}else{
return (0.0);}}
double DISTANCE(i,j)
int i;
int j;{
return (fabs(g1[i]-g1[j])+fabs(g2[i]-g2[j])+fabs(g3[i]-g3[j])+fabs(g4[i]-g4[j])+fabs(g5[i]-g5[j]));}double COUNT(epsilon,n)double epsilon;
int n;{
int i,j;double sum;for(i=0,sum=0.0; i<n; i++){for(j=0; j<n; j++){sum = sum + STEP(epsilon-DISTANCE(i,j));}}return ((double)sum/(n*n));}
ここは 81 ページ目カオス・フラクタル 2011 #9

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

main(){FILE *fpr;
int i,k;double x,k1x,k2x,k3x,k4x,kx,y,k1y,k2y,k3y,k4y,ky,z,k1z,k2z,k3z,k4z,kz;
double arg1,arg2,epsilon;
if((fpr = fopen("m5.dat", "wt")) !=NULL){
for(i = 0,x=0.01,y=0.01,z=0.01; i < imax; i++){k1x = h*func1(x,y,z);k1y = h*func2(x,y,z);
k1z = h*func3(x,y,z);k2x = h*func1(x+0.5*k1x,y+0.5*k1y,z+0.5*k1z);
k2y = h*func2(x+0.5*k1x,y+0.5*k1y,z+0.5*k1z);
k2z = h*func3(x+0.5*k1x,y+0.5*k1y,z+0.5*k1z);
k3x = h*func1(x+0.5*k2x,y+0.5*k2y,z+0.5*k2z);
k3y = h*func2(x+0.5*k2x,y+0.5*k2y,z+0.5*k2z);
k3z = h*func3(x+0.5*k2x,y+0.5*k2y,z+0.5*k2z);
k4x = h*func1(x+k3x,y+k3y,z+k3z);
k4y = h*func2(x+k3x,y+k3y,z+k3z);
k4z = h*func3(x+k3x,y+k3y,z+k3z);
kx = (k1x+2.0*k2x+2.0*k3x+k4x)/6.0;
ky = (k1y+2.0*k2y+2.0*k3y+k4y)/6.0;
kz = (k1z+2.0*k2z+2.0*k3z+k4z)/6.0;
x = x + kx;
xx[i]=x;y = y + ky;z = z + kz;}
for(k=0; k< N; k++){g1[k] = xx[k];g2[k] = xx[k+200];
g3[k] = xx[k+400];
g4[k] = xx[k+600];
g5[k] = xx[k+800];}
for(epsilon= 1.0; epsilon < 100; epsilon=epsilon+0.1){
arg1 = COUNT(epsilon,N);
arg2 = epsilon;
fprintf(fpr,"%lf %lf\n", log(arg2), log(arg1));}}
fclose(fpr);}
ここは 82 ページ目カオス・フラクタル 2011 #912.4

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

相関次元の一般化

前節で学んだ相関次元は単一のパラメータ q を導入することで次のように一般化することができる.

Dq∑n()log i=1 pqi1lim
q − 1 →0log =(198)

ここで, n() は半径  の球で対象 (軌道) を覆うために必要な球の個数であったから, 前節で調べた相関時限の場合には N = n() であることに注意しよう. そこで, 上式で q = 2 とおいてみると

D2=limlog∑N

2i=1 pi(199)log →0となり, 相関次元が復活する. また, q = 0 とおけば

D0=− lim

log

→0∑n()i=1 1log = − lim

→0log n()log (200)となり, 容量次元に一致する. また, q = 1 の場合も別な次元が定義される. ここで q → 1 の極限は0/0 を与え, 厄介に見えるが, 分母/分子を q で微分した後に q → 1 極限をとっても良いことに気が∑n()
つけば (いわゆる「ロピタルの定理」), 確率の規格化:
i=1 pi = 1 に注意して
} n(){∑n()n() q∑∑ qplogp∂iii=1pi log pilim=(201)log logpi  = lim∑n()q→1 ∂qq→1pq
i=1i=1i=1lim

q→1∂(q − 1)∂qi= 1(202)であるから∑n()lim Dq = D1q→1=lim

→0i=1 pi log pi(203)log が得られる. これを情報理論の平均情報量 (エントロピー) との形式的類似性から情報次元と呼ぶ.
従って, Dq は相関次元 D2 , 容量次元 D0 , および情報次元 D1 をその特別な場合として含む「一般化された次元」を与えている. 既に述べたように, これらの次元は場合によっては非整数次元となり得て, 再来週から学習していくことになる「フラクタル図形」を特徴つける際にも用いられる.

レポート課題 9
ローレンツ方程式 ( a = 10, b = 8/3, c = 30) に対し, m = 5, τ = 2 と選んだ場合の相関次元を計算せよ.

[カオス編] はこれでおしまいです. 次回 (7/5) から [フラクタル編] に移ります.

ここは 83 ページ目