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

平成 23 年 4 月 19 日

目次
4

非線形写像

4.1

ロジスティック写像から得られる数列のグラフを用いた求め方 . . . . . . . . . . .

4.2
4.3

初期値依存性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.4
4.5

非線形写像の周期について . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

数値計算上の丸め誤差/精度 (桁) 落ちによる影響 . . . . . . . . . . . . . . . . . . .

xn+1 = 4xn (1 − xn ) の周期解とその計算 . . . . . . . . . . . . . . . . . . . . . . .

131313141618

課題 1 の解答例
前回のこの講義ではロジスティック方程式の解析解, 2 つのタイプの差分化方程式を学んだ. 線形化20µ = 1, N0 = 10, λ = 0.1µ = 1, N0 = 10, λ = 0.1µ = 0.1, N0 = 20, λ = 0.5µ = 0.1, N0 = 20, λ = 0.5µ = 0.1, N0 = 2, λ = 0.5µ = 0.1, N0 = 2, λ = 0.1181614N12108642002t46810

図 5: ロジスティック方程式の解析解と線形化されたロジスティック方程式の差分方程式の比較.
されたロジスティック方程式の差分方程式とロジスティック方程式の解析解との間には各パラメー
タ (λ, µ, N0 ) の選び方によっては図 5 に示したような差異が見られるが, この誤差は評価すること11カオス・フラクタル 2011 #2

担当: 大学院情報科学研究科 井上 純一ができ, 必要とあらば改善することができる. それに関しては後日, この講義の「常微分方程式の数値計算」のところで系統的に学習することになる.
一方, 非線形のロジスティック方程式を差分化したものはこれらの結果と異なるのか否かを見ておくことは必要であろう. 非線形のロジスティック方程式を差分化したものにおいて, 例えば

λ = 0.1, µ = 1, N0 = 10 とおくと
Nt

2= (1 + 0.1)Nt−1 − 0.1Nt−1(33)であるから, はじめの数項を計算することにより

N1=(1 + 0.1)N0 − 0.1N02 = (1 + 0.1) × 10 − 0.1 × 100 = 1(34)N2=(1 + 0.1) × 1 − 0.1 × 1 = 1(35)であるから, 以降では常に 1 である. これは, ロジスティック方程式の解析解, 線形化したロジスティック方程式の解析解とも異なる. また, 前回の講義ノートの最後で述べたようにそれぞれの方
程式で λ = a − 1, µ = a/(a − 1), a > 1 とおき, a = 3.9 と選んでプロットすると図 6 のようになる. この図より, 解析解は滑らかにある一定値に近づくが, 線形化されたロジスティック方程式の差10.80.60.4

N 0.20Logistic map
discrete-linear
con.-linear-0.2-0.4-0.602040t6080100

図 6: λ = a − 1, µ = a/(a − 1), a = 3.9 と選んだ場合のロジスティック方程式の解析解 (cont-linear) と線形化されたロジスティック方程式 (discrete-linear), 非線形ロジスティック方程式の差分方程式 (Logistic map) の比較.

分方程式の振る舞いはこれと異なり, 初期段階では振動しながら最終的にはある一定値に近づいている (これらは第 1 回の講義ノート (14)(19) 式の形からもある程度予想できる). この一定値はロジスティック方程式の解析解のそれとは異なっている. 一方, 非線形ロジスティック方程式の差分方
程式の振る舞いはこれら 2 つの方程式の振る舞いのどちらとも全く異なり, 極めて不規則な挙動を示すことがわかる. また, この振る舞いを実際に計算する前に第 1 回講義ノート (22) 式から予想することは難しい.
今回から数回に渡る講義では, 不規則な振る舞いをする非線形ロジスティック写像の差分方程式を扱い, その振る舞いを詳細に見て行くことにする.

ここは 12 ページ目

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

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

非線形写像

4

ここまでで述べてきたことから, 非線形なロジスティック方程式を差分化した方程式:

Nt=2(1 + λ)Nt−1 − µλNt−1(36)の振る舞いは, それ自体がとても興味深いので, 研究や学習の対象となる. そこで, ここからはこの場合 — ロジスティック方程式の非線形版の差分方程式 — に焦点を絞って話しを進めることにしよう. その際, (36) 式をそのまま扱うのではなく, Nt = (a/µ(1 − a))xt+1 なる変数変換を施すことで, 方程式を簡略化しておく. つまり, ここからは xn についての更新式:

xn+1= axn (1 − xn )(37)の振る舞いについて調べていくことにしょう. この非線形な方程式をロジスティック写像と呼ぶ14.1

ロジスティック写像から得られる数列のグラフを用いた求め方

方程式 (37) を初期値 x = x0 から逐次求めることは計算機を用いれば難しくない. しかし,

y = f (x) = ax(1 − x) と y = x のグラフを描き, (x, y) = (x0 , 0) からスタートし, (x0 , x0 ) (y = xと x = x0 の交点), (x0 , f (x0 )) (y = x0 と y = f (x) の交点), (f (x0 ), f (x0 )) (y = x と x = f (0 ) の交点),・・・のように順次これらの点を辿って行くプロセスを繰り返し, その順番に x 軸の値を読みとったものが, (37) 式から生成される数列の値 x0 , x1 , x2 , · · · となっている事実を用いれば, 視覚的にこの数列を求めていくことができる. この図より, a = 2.5, 2.9 の場合には最終的に一つの点に収束していく様子が見て取れる. また, a = 3.1 の場合には x の 2 つの値の間を交互に行き来きする周期運動を示すこともわかる. しかし, a の値を増やし, a = 3.9 にした場合, x の値の変化から明確な規則を読み取ることは難しい. この a = 3.9 の場合の写像の振る舞いは既に図 6 でも見たように非常に複雑である. よって, ここからはいくつかの観点からこの複雑な振る舞いについての詳細を見ていくことにする.

4.2

初期値依存性

例えば a = 3.9 の場合, 初期値を x0 = 0.1 選んだが, この際にみられる複雑な挙動が x0 の値をx0 = 0.1 から x0 = 0.10001 に変えた場合にどうなるのかを調べてみよう. このような微小量の違いからくる振る舞いの差異を前節で説明したグラフによる作図から見るのは難しいので, 計算機を用いて写像 (漸化式) を逐次解くことでこれを数値的に調べてみる. 結果を図 8 に載せよう. この図より, n = 7, 8 ステップまで両者は完全に重なっているが, それ以降, これら両者の振る舞いには食い違いが生じる. 従って, この結果より, a = 3.9 の場合のロジスティック写像から生成される数列は初期値の選び方に対して非常に鋭敏であることがわかる. また, 図 7 で a = 2.2, 2.9 の場合のロジスティック写像に対しては, これらの初期値鋭敏性はみられない (各自, 調べてみること).

1 一般には a = 4 の場合の方程式をロジスティック写像と呼ぶが, この講義では a > 1 の任意の a での方程式をこのように呼ぶことにする.

ここは 13 ページ目

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

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

0.6xn+1xn+1x
path

1.2

0.598

1

2.9 x(1-x)

1.4x
path

0.602
1.2

0.7

2.5 x(1-x)

0.604

1.4

0.65

10.596 0.598

0.6xn

0.602 0.604

xn+1

xn+1

0.596
0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.2

0.4

xn

0.6

0.8

0

1

1

0.6

0.8

0.6

00.65
xn

0.7

0.2

0.4

xn

0.6

0.8

1

0.2

0.4

xn

0.6

0.8

1

1
3.1 x(1-x)

3.9 x(1-x)

xx

path

0.8

0.8

xn+1

0.6

xn+1

0.6

path

0.4

0.4

0.2

0.2

000.20.4xn

0.60.8100

図 7: ロジスティック写像からの数列のグラフを用いた求め方. y = f (x) = ax(1 − x) と y = x のグラフを描き,(x, y) = (x0 , 0) からスタートし, (x0 , x0 ) (y = x と x = x0 の交点), (x0 , f (x0 )) (y = x0 と y = f (x) の交点),(f (x0 ), f (x0 )) (y = x と x = f (0 ) の交点),・・・のように順次このプロセスを繰り返し, 点を辿って行く. 左上から右下にかけて a = 2.5, 2.9, 3.1, 3.9 の場合. いずれの場合も初期値は x0 = 0.1 に選んである.

4.3

数値計算上の丸め誤差/精度 (桁) 落ちによる影響

ここまでの考察で a = 3.9 の場合のロジスティック写像は初期値依存性を有することがわかった.
このような非線形写像が初期値の選び方に鋭敏である場合には, 数学上の同値変形から得られる 2つの数式であっても, それらは数値計算においては異なる結果を与える可能性があることに注意しなければならない.
例えば, 数学上は全く同じである次の 2 つの方程式を考えよう.

xn+1= axn (1 − xn )

(38)

xn+1

axn − ax2n(39)=2 年生のときに受講したであろう「計算機プログラミング I・同演習」では, 計算量の観点から上記2 つの方程式の表現法に対する計算結果の差について学んだクラスもあったかもしれない. 例えば,金曜クラス (A グループ) では

http://chaosweb.complex.eng.hokudai.ac.jp/~j_inoue/PROG2011/ProgI2011_5.pdfのようなことを学習した.
すなわち, (38) 式において, 掛け算の回数は 1 回である一方, (39) 式においては 3 回であるから（足し算/引き算はぞれぞれ 1 回ずつ), 可能な限りこの回数の少ない表現を用いてプログラムを書いた方が計算量の観点からは得であるというわけである.

ここは 14 ページ目

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

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

10.8xn0.6

0.40.2x0 = 0.1x0 = 0.10001
002040n6080100

図 8: xn+1 = 3.9xn (1 − xn ) において, 初期値を x0 = 0.1 に選んだ場合と x0 = 0.10001 に選んだ場合の数列の振る舞い.

しかし, 方程式が非線形の場合, これには注意が必要である. そこで, その問題点を具体的にみるために, a = 3.9 = 39/10, x0 = 0.5 = 1/2 として, (38) から得られる数列のはじめの数行を「分数」を用いて正確に書き出してみると

x1=391521
858879801, x2 =, x3 =,···4016000
2560000000(40)であるが, 計算機上は分母を分子で割ることで得られる少数を小数点以下どこで打ち切るかによる誤差の影響が無視できなくなる. また, 上に示したはじめの数項から推測できるように, 分母分子を正確に表現するための数字の桁数は繰り返しの増加とともに指数的に増加するので, 計算機の変数の表現範囲によっては, その制限から来る精度落ちの効果も無視できなくなる2 .
ここで注目すべきは, これらの影響が出始める回数 n は表現 (39) 式を用いた場合と一般的に異なるという点である. 実際, (39) 式に対してはじめの数項を書き出してみると

x1=39
39 39 39 39 39, x2 =×−×× ,···40
10 40 10 40 40(41)となっていくので, 例えば, 上の x2 の第 1 項において分子・分母の大きくなるスピードは第 2 項と比べて遅く, 従って, 打ち切りの効果が顕著に出始める繰り返し回数 n は (38) 式に基づく繰り返しと違う. よって, この打ち切り効果が異なり始めてその後の数列の振る舞いは, ちょうど, その時刻で初期条件を異なるものと置き換えて漸化式を作り直した状況として考えることができるので, 前に見たように初期値の違う非線形写像に対する初期値鋭敏性がここでも効いてきて両者の差はどんどん大きくなっていくことになる. この「差」が生まれるのはどのくらい繰り返しを行った後なのか, は使用する計算機に依存する. 参考までに担当者の所有するノート型パソコンで計算した結果を図 9 に載せる. この図からわかるように, 初期条件が同じであっても, 両者の間の際立った食い違いはだいたい n ' 6700 あたりで出始める.

2

例えば, 32 ビット単精度 (C 言語では int で宣言された整数) で任意整数を表現する際, 計算機内部では任意の十∑31進整数 z を 2 進整数 (b31 , b30 , · · · , b0 ) として z =b 2i と表現するので, これで表現しうる最大の整数はi=0 i(b31 , b30 , · · · , b0 ) = (1, 1, · · · , 1), すなわち, z = 1 + 2 + 22 + 23 + · · · + 231 = 232 − 1 = 4294967295 である. 従って, この値を越える整数値が分子あるいは分母に現れる場合, その整数値は何らかの規則 (四捨五入, 切り捨て等) でそれ以下の整数に変換されてしまう.

ここは 15 ページ目

カオス・フラクタル 2011 #2x

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

11

0.8

0.8

0.6

0.6

x
0.4

0.4

0.2

0.2

Type-1
Type-2002040n60

Type-1Type-28006640100110.80.80.6

6660

6680

n6700

6720

6740

9980

10000

0.6xx0.4

0.4

0.2

0.2

Type-1
Type-206700

6720

6740

n6760

Type-1
Type-26780

6800

09900

9920

9940

n9960

図 9: 式 (38)(これが Type-1), 式 (39)(これが Type-2) の振る舞い. 初期条件は同じ x0 = 0.5. だいたい n ' 6700 あたりで食い違いが出始める.

4.4

非線形写像の周期について

a = 3.9 あるいは, a = 4.0 に選んだロジスティック写像は非常に複雑な不規則運動を見せるが,繰り返し回数を十分取ったときに得られる数列の中には, いくつかの周期解が含まれるかもしれない. そこで, ここからはその周期について詳しくみていくことにする.
一般に非線形写像: xn+1 = f (xn ) に対し, ある x に写像を k 回繰り返し施したものを f k (x) で表すとすれば3 , f k (x) = x が成り立つとき, 写像 f は周期 k をもつと言う. 例えば, f 6 (x) = x ならば写像 f は周期 6 をもち, 関係式

f 6 (x) = f 3 (f 3 (x)) = f 3 (x) = x(42)から, f が 6 周期をもつならば, 3 周期をもつことになる. 従って, 一般に p を n の約数とすると,

f n (x) = x を満たす x の中には f p (x) = x を満たす x も含まれている. これを具体的にみるために, f (x) = 4x(1 − x) として, f k (x) = x を満たす x の値を調べてみよう. これは前節までに見てきたロジスティック写像の a = 4 の場合に相当する.
多少面倒ではあるが, 具体的に計算してみると

f (x) =4x(1 − x)

2f (x) =4{4x(1 − x)}(1 − {4x(1 − x)}) = 16x − 80x + 128x − 64x(44)

f 3 (x) =···(45)

(43)
234

3 「f の k 乗」ではないことに注意.

ここは 16 ページ目

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

担当: 大学院情報科学研究科 井上 純一であるから, 図 10 に示すように, 非線形写像 xn+1 = f (xn ) = 4xn (1 − xn ) の 1-4 周期解はそれぞ
れ y = f, · · · , f 4 と y = x の交点の x 座標の値である. この図からわかるように, 例えば, y = f 2 と110.80.80.6

0.6

0.4

0.4

0.2

0.2

000.2

0.4

x

0.6

0.8

10110.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

000.2

0.4

x0.6

0.8

1000.20.4x 0.6

0.8

100.2

0.4x 0.6

0.81

図 10: y = f 1 , · · · , f 4 と y = x. f (x) = 4x(1 − x) である. 左上から右下へ f 2 , f 4 , f 1 , f 3 .

y = x の交点は y = f と y = x の交点である 1 周期点 2 個を除いた 2 つの点からなり (2 は 1 を約数として持つ).
ところで, f 2 (x) = x を満たす 2 周期解 x を具体的に求めようとするのであれば, x に関する代数方程式:

15x − 80x2 + 128x3 − 64x4=0(46)を解かなければならない. この中には x = 0, 3/4 の 1 周期解も含まれているから, この左辺はx(x − 3/4) で割り切るはずなので, その商である x に関する 2 次方程式を解いて残りの 2 つの解を求めればよい. しかし, 一般的に n 周期解で n が大きくなると, この手の代数方程式の解を求めることは極めて厄介になる. ただし, a = 4 のロジスティック写像：xn+1 = f (xn ) = 4xn (1 − xn ) に限っては, 三角関数の倍角公式を使うことで, 周期解を首尾よく求めることができる. 以下にそれを見ておこう.

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

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

xn+1 = 4xn (1 − xn ) の周期解とその計算三角関数の倍角公式を用いるために, xn+1 = 4xn (1 − xn ) において xn = sin2 θn なる変数変換4を行うとsin2 θn+1= 4 sin2 θn (1 − sin2 θn ) = 4 sin2 θn cos2 θn = sin2 2θn(47)であるから, 結局, xn+1 = 4xn (1 − xn ) はこの変換によってsin2 θn+1

=sin2 2θn(48)と書き換えることができる. 従って, sin2 θn+2 が求めたいのであれば sin2 θn+2 = sin2 2θn+1 に対して上記を繰り返し使うことでsin2 2θn+1= 4 sin2 θn+1 cos2 θn+1 = 4 sin2 θn+1 (1 − sin2 θn+1 ) = 4 sin2 2θn cos2 2θn = sin2 4θn(49)であるから

sin2 θn+2=sin2 4θn(50)が得られる. 従って, 今欲しいのは p 周期解なので, これを θn+p に対して一般化すると

sin2 θn+p= sin2 2p θn(51)となる5 . 従って, これを満たすためには, m を自然数としてθn+p± 2p θn + mπ=(52)が成り立つべきである. 一方, p 周期である条件 (xn = sin2 θn = xn+p = sin2 θn+p ) から, やはり lを自然数としてθn+p = ± θn + lπ(53)が成り立つ. (52)(53) 式より, L ≡ l − m として

θn=Lπ(54)2p ± 1が得られる. 自然数 L には任意性があるが, ここでは簡単のため L = 1 と置くことにする.
まずは, 1 周期解 (x = 0, 3/4) が再現されるか確認してみると, p = 1 とおいて, θn = π, π/3 であるから, xn+p = sin2 θn+p = sin2 2p θn に注意すれば, それぞれから xn = sin2 π = 0, xn = sin2 π/3 =√( 3/2)2 = 3/4 となり, OK である. また, 2 周期解については, p = 2 として θn = π/3, π/5 が得ら√れるが, θn = π/3 のときには, xn = sin2 π/3 = ( 3/2)2 = 3/4 であり, これは 2 周期の中にその約数である 1 周期が含まれていることを表している. θn = π/5 のときには xn = sin2 π/5 = 0.345,また, (51) 式から xn+p = sin2 θn+p = sin2 2p θn であるから, xn+1 = sin2 2π/5 = 0.904, xn+3 =sin2 4π/5 = sin2 (π − π/5) = xn . よって, 2 周期解は 0.345 と 0.904 であり, これは図 10 からも見てとれる.
4 この変数変換を用いた議論は次回以降も使うので覚えておくと良い.
5 必要であれば帰納法で証明せよ.

ここは 18 ページ目

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

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

ここで見た周期解は安定であるわけではなく, 不安定化する. その安定性条件に関する詳細は次回以降に見ていくことにする.

レポート課題 2
複雑な振る舞いを見せるのは何もロジスティック写像だけではない. 例えば次のようなテント写像:

xn+1= 1 − |1 − 2xn |(55)

もロジスティック写像同様に複雑な挙動を示す. この写像のはじめの 100 ステップの値を初期値

x0 = 0.1 に対して計算し, 横軸に繰り返し回数 n, 縦軸に xn の値をプロットせよ6 . この結果は次回の講義内容に関係するので, 各自が必ず計算, グラフ化し, 結果を目で見てみること.

連絡: 次週 4/26 は休講とします.

6

どうもうまく行かない場合には 2 = 1.9999 として, xn+1 = 1 − |1 − 1.9999xn | を用いよ.

ここは 19 ページ目