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

平成 23 年 5 月 24 日

目次
8

分岐現象とカオス
8.1 ロジスティック写像 La (x) 再考 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4242

固定点とその安定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8.2

8.1.1
8.1.2

2 周期解とその安定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43448.1.3
8.1.4
8.1.52n 周期解とその安定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
カオスが出現する分岐パラメータ a 値の粗い見積もり . . . . . . . . . . . .
窓, 自己相似構造 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

464647

参考: ニュートン法の復習 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

課題 4 の解答例
まずは前回学んだリアプノフ指数の計算例を見ておこう. テント写像はパラメータ a を導入することで次のように一般化することができる.

Ta (x)= a(1 − |1 − 2x|), xn+1 = Ta (xn )

(103)
0

リアプノフ指数の定義から, (103) 式に基づいて生成される軌道の各点での |Ta (x)| が求まればリアプノフ指数が計算できる. これは明らかに x の値に依らず 2a であるから数値計算するまでもなく
01∑n log 2a= log 2alog |Ta (xi )| = limn→∞ nn→∞ni=0n−1λa=lim(104)である. ということは, a を 1/2 ≤ a ≤ 1 の範囲で変化させると必ずリアプノフ指数は正の値になり, パラメータ a で一般化されたテント写像から生成される軌道はカオスとなる. また, a を変化しつつ, 写像からの軌道の取りうる値を重ね打ちすると興味深い現象が見られる. その結果をリアプノフ指数の a 依存性とともに図 22 に載せる.
参考までにロジスティック写像 La6=4 (x) のリアプノフ指数を計算するためのプログラムの一部を以下に載せておこう1 .

/* 写像が陽に与えられた場合のリアプノフ指数の数値計算 */for(a=1.00; a <= 4.00; a = a + 0.0001){
1 ここではファイルのオープン/クローズを使わないが, ファイル名を指定してリダイレクションすればよい.

41

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

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

1generalized tent map
log(2*x)0.9
0.8
0.7x0.6
0.5
0.4
0.3
0.2
0.1
00.5

0.6

0.7

a0.80.91

図 22: 一般化されたテント写像の分岐図と対応するリアプノフ指数.

for(i = 0,sum=0.0,x=x0; i <= N; i++){// x = a*(1.0-fabs(1-1.999*x)); /* テント写像ではこちらを使う */x= a*x*(1-x);
sum = sum + log(a*fabs(1.0-2.0*x)); /* 写像の出力結果 x に対して和を計算 */}printf("%lf %lf\n",a,(double)sum/N);}

分岐現象とカオス8

今回の講義では今まで学んできたロジスティック写像 La (x) に対し, パラメータ a の値を変えていくことで, 写像からの軌道が周期運動からカオスへと段階的に変化する様子を「分岐現象」として捉え, その分岐が生じるたびに今まで安定であった周期解が不安定化し, さらに周期の長い解が安定解として現れる様子を見ていく. a の変化による周期解の不安定化と, さらなる周期解の出現を繰り返し, ある臨界パラメータ a 以上でカオスが現れることを学習する. また分岐図には「窓」や「自己相似構造」などの特異な性質が現れることも同時に見ていく.

8.1

ロジスティック写像 La (x) 再考

まずは図 23 にロジスティック写像のパラメータ a に対する軌道の分岐図を載せよう. この図は横軸にパラメータ a を縦軸に写像: xn+1 = axn (1 − xn ) に対し, xn の取りうる値を重ねうちしたものである2 . 先に載せたリアプノフ指数の計算とほぼ同じループ処理を行えば良いのだが, 念のため下記に図 23 を計算する際に用いたプログラムの一部を載せておく.

2 細かな計算技術上の話になるが, 固定点に収束する場合でさえも, その固定点へ 1 ステップで収束するとは限らない.

多くの場合, このような固定点 (あるいは周期解) へ収束するまでには有限回のステップを必要とする. このような収束までに経過する x の点までプロットするとうまくないので, はじめの数ステップはかさ値打ちせずに捨てるなどの工夫をすると良い.

ここは 42 ページ目

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

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

10.90.8
0.7
0.60.50.4
0.30.20.1011.522.5a33.54

図 23: ロジスティック写像のパラメータ a に対する軌道の分岐図.

/* 分岐図計算のプログラム */for(a=1.00; a <= 4.00; a = a + 0.0001){
for(i = 0,sum=0.0,x=x0; i <= N; i++){// x = a*(1.0-fabs(1-1.999*x)); /* テント写像の場合はこちらを使う */x= a*x*(1-x);
if(i >= Nskip){
printf("%lf %lf",a,x);
}}
}

8.1.1

固定点とその安定性

この図より, パラメータ a を a = 1 から徐々に大きくしていくと a = 3 付近で固定点は 2 周期軌
道に分岐する. そこで, まずは固定点の位置 x∗ を a の関数として求めるために, La (x) = x を x について解いてみると, x∗ = 0, および, x∗ = (a − 1)/a が得られる. これを図 23 に重ねうちしたものを図 24 に載せる. この図より, a = 3 近傍まで x∗ = (a − 1)/a と写像 xn+1 = La (xn ) から得られる軌道のプロットは完全に一致するが, a > 3 では実際に写像から得られる軌道は x∗ = (a − 1)/aではなく, その両側に現れる「2 周期軌道」となっている. 従って, 固定点 x∗ = (a − 1)/a は a > 3においても解としては存在するが, それは安定ではない.
そこで, リアプノフ指数の導出の議論で見たように, |{dLa (x)/dx|x=x∗ }| が 1 をこえる条件を調べてみると, 直ちに{dLa (x)dx}{=x=x∗2(a − 1)a 1−a}= |2 − a| > 1(105)が得られ, 従って, これは (1 − a)(3 − a) > 1 すなわち, a > 3 である. よって, 固定点 x∗ = (a − 1)/aは a > 3 で不安定化する (1 < a < 3 で安定である).
ここは 43 ページ目

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

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

10.8map1-(1/a)
0.6x
0.4

0.2

0

1

1.5

2

a2.5

33.5

図 24: 写像から得られる分岐の様子と固定点を a の関数として求めたもの x∗ = (a − 1)/a を重ねうちしてプロットした.
ところで, La (x) = x のもう一つの解 x∗ = 0 は実現されないのか, 気になるかもしれない. 実際,この解は常に存在するのであるが, これもやはり安定に実現されるための条件がある. そこで, 上記のように La (x) の一階微分の x∗ = 0 での絶対値が 1 をこえる条件 (不安定化条件) を出してみると |{dLa (x)/dx|x=0 }| = a > 1 であるから, a > 1 の領域で x∗ = 0 は決して実現されない. 逆に,a ≤ 1 であれば x∗ = 0 は安定な固定点である (各自が写像の反復数値計算で確かめてみよ).

8.1.22 周期解とその安定性

図 23 の数値計算結果によると a > 3 では固定点の代わりに 2 周期解が実現されるわけであるから, この 2 周期解を求めてみる. これは L2a (x) = La (La (x)) = x を解けばよい. x, (−ax + a − 1)を因数に持つことに注意して (x∗ = 0, (a − 1)/a は常に解として存在する)

L2a (x) − x =x(−ax + a − 1){a2 x2 − a(a + 1)x + a + 1} = 0(106)すなわち, a2 x2 − a(a + 1)x + a + 1 = 0 の実解である√√
a + 1 + (a + 1)(a − 3)
a + 1 − (a + 1)(a − 3)(+)(−)
x∗2 =, x∗2 =2a2a(107)(±)の 2 つが 2 周期解の「枝」であり, それぞれ上記のような a の関数である. 従って, x∗2 が安定か(±)

どうかは L2a (x) の x∗2 における微分係数が 1 よりも大きいかどうかを調べればよい3 . 前回見た陰関数の微分を思い出すと

dL2a (x)dx=dL2a (x) dLa (x)·= a(1 − 2La (x)) · a(1 − 2x)
dLa (x)dx(108)(±)

であるから, この右辺を x∗2 で評価すればよい.

(+)

これは面倒で骨の折れる作業であるが, 以下で x∗2 について地道に計算してみると, まず(+)

La (x∗2 ) =
3(+)

(+)

ax∗2 (1 − x∗2 )2 周期解の安定性を調べるのであるから, La (x) の微分係数ではなく L2a (x) の微分係数を見る必要があることに注意.

ここは 44 ページ目

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

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

1map(a+1+sqrt((a+1)*(a-3)))/(2*a)(a+1-sqrt((a+1)*(a-3)))/(2*a)

0.8

0.6

x0.4

0.201

1.522.5a33.54(−)
図 25: 写像から得られる分岐の様子と 2 周期点を a の関数として求めたもの x(+)
∗2 , x∗2 を重ねうちしてプロットした.

√√(a + 1)(a − 3) a − 1 − (a + 1)(a − 3)=·2a√ 2
a + 1 − (a + 1)(a − 3)(−)
= x∗2=2a√(+)
が得られ, 従って, a(1 − 2La (x∗2 )) = (a + 1)(a − 3) − 1 である. 一方√(+)
a(1 − 2x∗2 ) = −1 − (a + 1)(a − 3)
a+1+(109)(110)であるから, 結局

dL2a (x)dx√√= −( (a + 1)(a − 3) − 1)( (a + 1)(a − 3) + 1) = −(a2 − 2a − 4) (111)(+)

x=x∗2となる. 従って, |a2 − 2a − 4| < 1 が安定条件であるが, 両辺を自乗して, 右辺の 1 を左辺に移し√√てまとめると (a + 1)(a + 6 − 1)(a − 3)(a − 1 − 6) < 0 が得られる. よって, a − 3 > 0 かつ√√(+)
a − 1 − 6 < 0, すなわち, 3 < a < 1 + 6 で x∗2 が安定に存在する4 .
(−)

(−)

(+)

次に x∗2 の安定性であるが, L2a (x∗2 ) = x∗2 であるから√(−)
a(1 − 2L2a (x∗2 )) = −1 − (a + 1)(a − 3)(112)であり, また(−)

a(1 − 2x∗2 ) =−1 +√(a + 1)(a − 3)

(113)であるから, (108) 式より

dL2a (x)dx= −(a2 − 2a − 4) =(−)

x∗2(−)(+)

dL2a (x)
dx

(114)
(+)

x∗2が得られる. 従って, x∗2 の安定領域は x∗2 の安定領域と同じ 3 < a < 1 +√

6 であることがわかる. これは図 23 からも見て取れる.
4

はじめの 2 項の積: (a + 1)(a +√6 − 1) はこの領域で常に正なので問題にはならないことに注意.

ここは 45 ページ目

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

8.1.3

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

2n 周期解とその安定性√6 を越えると不安定化し, a > 1 + 6 での安定解は 2 周期解から 4 周期解に移る. 従って, 前節までの議論をここで繰り返すためには L4a (x) = x を(±)
満たす x = 0, (a − 1)/a, x∗2 以外の解を求め, その解に対する微分係数 dL4a (x)/dx の絶対値を計(±)

図 23 を見る限りにおいて, x∗2 も a = 1 +√算し, それと 1 と比較すればよい.
図 23 の分岐図から, その後, a の増加とともに周期は倍々に増えていくことがわかるので, その解の安定性の議論も上記と同じ手続きをふめばよい. しかし, 系統的な手続きがわかっていることと, それが実際に実行可能であるか否かは別問題である. 例えば, L4a (x) = x を具体的に書き下してみるとa(a(a(ax(1 − x))(1 − (ax(1 − x))))(1 − (a(ax(1 − x))(1 − (ax(1 − x))))))
×(1 − a(a(ax(1 − x))(1 − (ax(1 − x))))(1 − (a(ax(1 − x))(1 − (ax(1 − x)))))) − x = 0(115)(±)となり, x に関する 16 次方程式となる. 上記方程式は x = 0, (a − 1)/a, x∗2 を解に持つので,(+)(−)

x(x − (a − 1)/a)(x − x∗2 )(x − x∗2 ) を因数に持つことがわかるから, この多項式で (115) 式左辺を割れば, その商である方程式の次数は 12 次まで減少する. しかし, これでもその実解を解析的に求めることは不可能である.
このような場合には, 2 年生のときに履修した「計算機プログラミング I」の講義で学んだであろう「ニュートン法」などの数値計算技法を用いて多項式のゼロ点を探すことになる (⇒ 今週のレポート課題).

8.1.4

カオスが出現する分岐パラメータ a 値の粗い見積もり

前節で述べたように, 周期はパラメータ a の増加とともに倍々で増加してゆき, ある極限でカオスが現れることが予想される. この分岐の詳細を正確に議論することは前述の理由で難しいのであるが, かなり粗い見積もりならばカオスが出現する a の近似値を評価することができる. 以下でそれを述べよう.
今, 我々が正確にわかっているのは x = 0 から x = (a − 1)/a への分岐点: a0 = 1, x = (a − 1)/a(±)(±)

から, x = x∗2 への分岐点: a1 = 3, x = x∗2 から 4 周期解 (具体的には求まっていないが) への分√岐点: a2 = 1 + 6 である. そこで, a2 − a1 , a1 − a0 の比, つまり, 「階差の比」を計算してみると

a2 − a1a1 − a0=√√
1+ 6−36−2== 0.225 ≡ r
3−12(116)となる. 次の階差比: (a3 − a2 )/(a2 − a1 ) を計算するためには a3 の具体的な値が必要であるが,我々は (さしあたり) それを求めることができない. そこで, これ以降, 任意の階差比: (an+2 −an+1 )/(an+1 − an ) についても (116) 式の公比 r が有効であると仮定する. すなわち, 階差は次の等比級数に従うとする.

an+2 − an+1=r(an+1 − an )(117)

これを具体的に書き下してみると

a2 − a1=r(a1 − a0 )ここは 46 ページ目カオス・フラクタル 2011 #5担当: 大学院情報科学研究科 井上 純一

a3 − a2=r(a2 − a1 )··· ··· ···a∞ − a∞−1=r(a∞−1 − a∞−2 )となるから, これらを辺々足していくと

a∞ − a1= r(a∞−1 − a0 )(118)が得られる. そこで a∞ = a∞−1 と置き, a∞ について解くとa∞=a1 − ra0
3 − 0.225 × 1== 3.58
1−r
1 − 0.225(119)となる. 数値計算による正確な値は a∞ = 3.569 であることがわかっているので5 , これと比べてもかなり良い近似値になっていることがわかる. この結果と第 4 回講義ノート図 18(右) のリアプノ
フ指数が正となる領域の対応を考えてみるとよいであろう.

8.1.5

窓, 自己相似構造

a > a∞ のパラメータ領域に対する分岐図を一部拡大してみると図 26 が得られる. この図からわかるように, ところどころに軌道が密でない周期解を持つ「窓」(4 周期窓や 3 周期窓) が見られる. また, この中で 3 周期窓の領域を越えたところを着目すると, 図 23 の縮小版のような分岐図が見られ, 自己相似な「入れ子構造」を見ることができる. このような微細な構造の詳細については,時間があれば, この講義の後半 (フラクタル編が一通り終わった後) で見ていくか, 時間が無ければ,後期「情報工学演習 II」でとりあげたい.

8.2

参考: ニュートン法の復習

f (x) = a2 x2 − a(a + 1)x + a + 1 = 0 の手の方程式 (x の多項式 = 0) の数値解を求めたい場合,ニュートン法と呼ばれる手法を用いるとよい.
この方法は方程式 f (x) = 0 の解の候補 x = x0 で曲線 y = f (x) に接する接線の方程式:
y=0f (x0 ) + f (x0 )(x − x0 )(120)の x 軸との交点 :x1≡ x0 −f (x0 )f 0 (x0 )(121)を解候補 x0 の改良版とするもので, これを推し進めて n-番目の解候補と n + 1-番目の解候補の間に成り立つ漸化式

xn+1= xn −f (xn )f 0 (xn )(122)を反復的に解き, この収束点を f (x) = 0 の解とするものである.
※注:「2 分法」というやり方でも解を求めることができる. 詳しくは5

この正確な値を求めるための, やや込み入った方法は後期開講の「情報工学演習 II」で扱う (多分).

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

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

10.8

0.6

x0.4

0.2

03.6

3.65

3.7

3.75

3.8a3.85

3.9

3.95

40.1
3.82

3.825

3.83

3.835

3.84

a 3.845

3.85

3.855

3.86

1
0.90.80.7

x0.6
0.5
0.4
0.30.2

図 26: ロジスティック写像のパラメータ a に対する軌道の分岐図. 図 23 の一部を拡大したもの.

http://chaosweb.complex.eng.hokudai.ac.jp/~j_inoue/PROG2010/ProgI2010_5.pdf
を参照してみるとよい.
レポート課題 52 次方程式:
f (x) ≡ a2 x2 − a(a + 1)x + a + 1=0をニュートン法により求めるプログラムを作成せよ. このとき

• x の初期値 x0 を [0, 1] の一様乱数により与え, 初期値をいくつか変えることで全ての実解が自動的に求まるようなプログラムを組むとなお良い.
√• a の値を 3 < a < 1 + 6 の範囲内で変化させると, 解の値も変化するが, この解を a の関数として求めるようなプログラムにするとなお良い.
• 求めた解 x に対し, |dL2a (x)/dx| − 1 の値を計算し, これがゼロ以上か否かを判定する機能を持たせるとなお良い.
上記 3 つの機能は「オプション」. 初期値を手で与えて解を求めるプログラムが書ければよい.

ここは 48 ページ目